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ABSTRACT 


In  this  research  a  multiple  finite  source  queueing  model 
with  a  single  server  and  fixed  priority  service  discipline  is 
investigated.  Interarrival  times  have  exponential  density  func¬ 
tions  and  service  times  have  arbitrary  distribution  functions . 

A  solution  procedure,  containi'-.'  a  set  of  numerical  algorithms, 
is  developed  to  obtain  server  utilization  and  the  mean  values 
of  the  waiting  times  per  customer  of  different  classes  and  the 
number  of  waiting  customers  of  different  classes. 

The  imbedded  Markov  chain  technique  is  used  as  the  basic 

method  of  analysis  of  the  model.  First,  the  expressions  for 

the  utilization,  p,  and  the  mean  values  of  the  system  performance 

measures  for  different  classes  are  obtained  in  terms  of  p.'s,  the 

1 

proportion  of  time  the  server  is  busy  with  the  customers  of  the 
corresponding  classes,  using  an  extension  of  Little's  formula. 
Based  on  regenerative  theory,  these  p^'s  are  related  to  ECB^'s, 
(i  =  1,2,..., K),  which  are  the  expected  lengths  of  the  busy 
period  during  which  the  server  is  busy  with  class  i  customers 
and  E(I),  the  expected  length  of  the  idle  period.  The  E(B^)'s 
are  then  expressed  as  functions  of  the  steady  state  probabilities 
of  the  Markov  chain  obtained  by  imbedding  the  process  at  depar¬ 
ture  epochs. 

After  the  elements  of  the  transition  probability  matrix  of 
the  Markov  Chain,  P,  are  generated,  obtaining  the  steady  state 


probabilities  involves  inversion  of  a  large  matrix.  Numerical 
techniques  are  developed  for  inversion,  first  by  reordering  the 
states  of  P  to  have  a  suitable  structure  for  the  matrix  to  be 
inverted,  and  then  modifying  an  available  inversion  procedure 
to  suit  this  structure.  Using  the  expressions  developed  for  the 
inverse,  recursive  relations  are  derived  for  finding  the  values 
of  steady  state  departure  probabilities. 

The  set  of  algorithms  developed  are  extended  to  find  time 
average  marginal  and  joint  probabilities  of  finding  a  certain 
number  of  waiting  customers  of  different  classes  at  the  facility 
using  Level  Crossing  Analysis.  The  possibility  of  using  the 
algorithms  fox  mixed  class  models  consisting  of  finite  sources 
for  some  classes  and  infinite  sources  for  the  other  classes  is 
also  discussed. 
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CHAPTER  1 


INTRODUCTION  AND  REVIEW  OF  THE  LITERATURE 

1.1  INTRODUCTION 

Queueing  theory  serves  as  a  useful  mathematical  tool  to 
analyze  the  effects  of  various  queueing  phenomena  to  be  con¬ 
sidered  in  the  modeling  and  analysis  of  many  systems.  A  mathe¬ 
matical  study  of  any  system  using  queueing  models  needs  speci¬ 
fication  about  the  capacity  of  the  source  from  which  customers 
are  generated,  the  distributions  of  the  inter-arrival  and  ser¬ 
vice  times  of  the  customers,  the  number  and  arrangement  of 
servers  and  service  stations  at  the  service  facility,  and  the 
service  discipline  based  on  which  the  customers  are  selected 
for  service.  When  a  mathematical  model  of  a  system  Is  construc¬ 
ted,  the  underlying  motivation  is  usually  to  evaluate  some 
measure  of  performance.  In  the  case  of  queueing  systems, 
usually  the  performance  measures  which  are  of  interest  are 
the  waiting  times,  the  number  of  customers  in  the  queue  and 
at  the  service  facility,  and  the  utilization  of  the  server. 

The  mean,  variance,  and  the  probability  distribution  of  these 
variables  are  studied. 

Queueing  models  in  which  customers  are  generated  from  at 
least  one  finite  input  source  to  which  they  return  after  re¬ 
ceiving  service  are  known  as  finite-source  queueing  models. 


Such  finite  source  models  were  initially  used  to  study  indus¬ 
trial  processes  in  which  one  operator  or  a  group  of  operators 
attend  a  finite  number  of  machines  which  may  break  down  from 
time  to  time.  Therefore,  in  queueing  theory  literature,  these 
models  are  known  as  the  machine  servicing  model  [FELL  66]  or 
the  machine  interference  model  [COXS  61].  Each  machine  is 
either  running,  requiring  repair  service,  or  waiting  as  a 
standby.  When  a  machine  breaks  down,  it  joins  the  queue  of 
machines  waiting  for  repair.  If  the  operator  is  free,  a 
machine  is  selected  from  the  queue  based  on  some  scheduling 
rule.  If  the  operator  is  busy,  the  machine  must  wait  in  the 
queue  until  it  is  selected  for  service.  After  a  machine  is 
serviced,  it  starts  working  and  after  a  passage  of  some  time 
interval,  which  may  be  stochastic,  fails  again  and  requires 
the  service  of  the  operator. 

Such  finite  source  models  are  also  models  of  many  time¬ 
sharing  computer  and  multi-access  communications  systems  in 
which  a  finite  number  of  users  or  computer  terminals  depend 
on  the  service  from  a  computer  system.  Each  user  sends  from 
the  terminal  a  request  for  processing  to  the  computer  and  the 
system  keeps  this  request  in  the  queue.  When  the  particular 
request  is  selected  for  processing  according  to  the  scheduling 
rule  used  for  the  processor,  the  program  associated  with  it 
is  executed.  After  the  execution  is  completed,  the  response 
or  the  output  is  fed  back  to  the  terminal  and  then  the  user 
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begins  to  generate  a  new  request  for"  the  computer. 


If,  in  such  models,  more  than  one  type  of  customer  ema¬ 
nates  from  one  or  more  input  sources,  then  the  question  of 
allocating  priorities  arises  and  is  of  practical  significance. 
Two  types  of  models  may  arise  in  such  cases  with  priorities 
[JAIS  67].  They  are  (i)  the  multiple  finite  source  priority 
models,  in  which  K  types  (K  >_  2)  of  customers  are  generated 
from  K  different  independent  sources,  and  (ii)  the  single 
finite  source  priority  models  in  which  K  types  (K  >_  2)  of  cus¬ 
tomers  are  generated  from  a  single  finite  source.  An  example 
of  the  multiple  finite  source  priority  model  is  the  multiaccess 
computer  system  which  is  used  by  different  independent  classes 
or  groups  of  finite  number  of  users.  The  priority  given  to 
a  user  when  being  selected  for  service  can  depend  on  the  class 
to  which  the  user  belongs.  Single  finite  source  models  are 
best  illustrated  by  the  situation  in  which  an  operator  looks 

after  a  set  of  N  number  of  machines,  each  of  which  can  fail 

because  of  any  of  K  types  of  failures,  with  certain  probabili¬ 
ties.  Based  on  the  type  of  failure,  some  priority  rule  is 
used  to  select  a  broken  down  machine  for  repair.  The  priority 

rule  used  in  a  queueing  model  can  be  of  different  types.  One 

of  the  most  commonly  used  priority  disciplines  assigns  dif¬ 
ferent  fixed  priorities  to  different  groups  or  types  of  cus¬ 
tomers  based  on  some  external  characteristics  and  always  gives 
preference  to  higher  priority  customers  over  those  with  lower 


priority.  This  implies  that  a  lower  priority  customer  is  taken 
for  service  only  if  there  are  no  higher  priority  customers 
present.  If  the  service  of  a  lower  priority  customer,  already 
being  serviced,  is  interrupted  before  completion  of  service  on 
the  arrival  of  a  higher  priority  customer,  then  this  priority 
discipline  is  called  preemptive  fixed  priority  discipline.  If 
the  service  of  a  customer  is  never  interrupted  before  completion, 
then  this  type  of  service  discipline  is  known  as  non-preemptive 
fixed  priority  discipline.  In  both  types  ,  the  customers  within 
the  same  class  are  selected  for  service  based  on  the  order  of 
their  arrival. 

Finite  source  queueing  models  are  difficult  to  analyze,  as 
compared  to  infinite  source  models,  because  the  arrival  rate 
of  the  customers  at  the  service  facility  is  not  constant,  but 
is  dependent  on  the  number  of  customers  already  at  the  service 
facility.  Only  limited  work  has  been  done  in  finite  source 
models  as  compared  to  infinite  source  models.  The  complexity 
of  the  analysis  increases  as  the  number  of  classes  increases. 
Systems  with  large  source  capacities  can  be  analyzed  approxi¬ 
mately  as  infinite  source  models  if  the  arrival  rates  of  the 
customers  at  the  facility  are  less  dependent  on  the  number  of 
customers  already  at  the  facility.  But  such  an  approximation  may 
not  be  valid  in  systems  in  which  the  source  capacities  are 
sufficiently  small  and  the  arrival  rates  of  the  customers  are 
dependent  on  the  number  of  customers  already  at  the  facility. 
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Multiple  finite  source  priority  models  are  comparatively 
easier  to  analyze  than  single  finite  source  priority  models. 
This  is  because  in  the  case  of  single  finite  source  priority 
models  the  arrival  of  any  class  customer  at  the  facility 
affects  the  arrival  rates  of  all  other  classes,  which  is  not 
the  case  with  multiple  finite  source  priority  models.  As  far 
as  the  priority  disciplines  are  concerned,  analysis  of  preemp¬ 
tive  priority  models  is  easier  compared  with  that  of  non-pre- 
emptive  priority  models,  because  the  lower  priority  customers 
do  not  have  any  effect  on  the  higher  priority  customers  in 
preemptive  discipline. 


1.2  OBJECTIVE  OF  THIS  RESEARCH 

In  this  research,  the  multiple  finite  source  priority 
model  with  a  single  server  and  non-preemptive  fixed  priority 
service  discipline  is  considered.  The  objective  is  to  develop 
numerical  algorithms  to  obtain  mean  waiting  times  of  each  class 
customer,  the  mean  numbers  of  each  class  customers  waiting  in 
the  queue  and  at  the  service  facility,  and  the  utilization  of 
the  server. 

The  schematic  diagram  of  the  model  considered  in  this  re¬ 
search  is  given  in  Figure  1.1.  There  are  K(>^  2)  classes  of 
customers,  originating  from  K  independent  finite  sources.  The 
capacity  of  source  i,  (i  =  1,2,...,K),  is  which  is  finite. 
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SCHEMATIC  DIAGRAM  OF  THE  MODEL 


Each  customer  of  class  i  calls  for  service  at  the  service  facil- 


.  i 


ity  after  spending  a  random  amount  of  time  at  the  source.  This 
time  is  exponentially  distributed  with  mean  1/A^.  There  is  a 
single  server  at  the  service  facility.  The  service  times  of 
class  i  customers  are  identically  and  independently  distributed 
random  variables,  with  an  arbitrary  distribution  function  FCs^) 
and  with  a  mean  1/p^,  where  is  the  mean  service  rate  of  class 
i  customers.  Non-preemptive  fixed  priority  service  discipline 
is  used  by  the  server  with  higher  priority  given  to  class  i 
customers  compared  with  class  j  customers  if  i  <  j .  Within  the 
same  class,  customers  are  selected  for  service  based  on  the 
order  of  their  arrival.  Extending  Kendall's  notation  and  the 
report  of  the  Queueing  Standardization  Conference  [MOOR  72], 
this  model  can  be  written  as  a  M/G/l/  /N^N^..., 
preemptive  discipline. 

There  are  several  methods  available  to  analyze  queueing 
models,  among  which  some  are  specially  suitable  for  the  model 
in  this  research.  These  are  described  in  the  next  section. 

1.3  METHODS  FOR  ANALYSIS  OF  THE  MODEL 

A  queueing  process  is  basically  a  stochastic  process  and  the 
method  of  analyzing  any  queueing  model  depends  on  the  type  of  the 
inherent  stochastic  process  of  the  model.  An  interesting  class 
of  stochastic  processes  is  the  Markov  process  in  which  the  pre¬ 
sent  state  of  the  system  is  sufficient  to  predict  the  future 
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without  a  knowledge  of  the  past  history  of  the  system.  In  case 
of  any  queueing  model  with  exponential  inter-arrival  and  service 
times,  the  number  of  customers  present  at  the  service  facility 
at  any  arbitrary  time  forms  a  Markov  process.  In  such  models, 
taking  the  state  of  the  system  as  given  by  the  number  of  custom¬ 
ers  present  at  the  service  facility,  the  equilibrium  or  the 
steady  state  time  average  probability  distribution  of  the  num¬ 
ber  of  customers  present  at  the  facility  can  be  obtained.  This 
is  done  by  writing  the  balance  equations,  which  are  known  as 
global  balance  equations,  for  each  state  based  on  the  principle 
that  the  probabilistic  flow  rate  into  a  state  must  equal  the  prob 
a bilistic  flow  rate  out  of  that  state  under  steady  state  condi¬ 
tions.  These  global  balance  equations  form  a  set  of  simultane¬ 
ous  equations  which  are  then  solved  to  find  the  steady  state 
time  average  probabilities.  In  some  cases,  it  is  possible  to 
derive  recursive  relations  between  the  probabilities  of  succes¬ 
sive  states  which  may  lead  to  a  general  equation  for  the  prob¬ 
abilities.  This  saves  time  required  for  solving  the  simultane¬ 
ous  equations. 

When  either  the  service  times  or  the  inter-arrival  times 
of  the  customers  in  a  model  are  not  exponential,  the  number  of 

customers  present  at  the  service  facility  at  any  arbitrary  time 
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is  no  longer  a  Markov  process.  There  are  three  commonly  used 
approaches  to  handle  such  cases  which  are  suitable  for  the 
model  considered  in  this  research  [KLEI  75].  These  are  dis- 


cussed  in  the  following  subsections  with  respect  to  the  expo¬ 
nential  arrival  time  and  arbitrary  service  time  distributions. 

1.3.1  The  Generalized  Method  of  Stages 

One  approach  to  solve  nonexponential  service  time  distribu¬ 
tions  is  the  method  of  stages  and  its  generalization.  Cox 
[COX  55]  showed  that  any  probability  density  function  having  a 
rational  Laplace  transform  can  be  represented  as  a  combination 
of  fictitious  exponential  stages  with  appropriate  mean  for  each 
stage.  The  advantage  in  such  a  representation  of  an  arbitrary 
service  time  distribution  is  that  by  including  the  stage  in 
which  the  customer  is  in  service,  in  the  description  of  the  state 
along  with  the  number  of  customers  at  the  facility,  the  global 
balance  equations  can  be  written. 

The  main  problem  is,  however,  that  the  dimension  of  the 
state  space  rises  rather  sharply  as  the  system  complexity  and 
the  number  of  stages  in  the  representation  of  the  service  time 
distribution  grow.  This  correspondingly  increases  the  number 
of  simultaneous  equations  to  be  solved,  though  it  is  easy  to  ob¬ 
tain  the  matrix  corresponding  to  these  equations.  Also,  this 
approach  is  obviously  restricted  to  service  time  distributions 
with  rational  Laplace  transforms. 
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1.3.2  Imbedded  Markov  Chain  Analysis 

Any  non-Markovian  process  can  be  studied  by  extracting  a 
set  of  points  at  which  the  Markovian  property  holds.  Such  points 
are  called  regeneration  points.  In  queueing  models  with  exponen¬ 
tial  inter-arrival  times  and  general  service  times,  the  epochs 
at  which  customers  depart  from  the  service  facility  constitute 
a  set  of  regeneration  points.  Therefore  the  number  of  customers 
present  at  the  service  facility  at  these  points  forms  a  Markov 
process.  This  method  is  called  the  imbedded  Markov  chain  analy¬ 
sis  because  it  involves  extracting  a  discrete-time  Markov  chain 
imbedded  in  the  continuous-time  process  and  this  technique  is  due 
to  Kendall  [KEND  50]. 

Using  the  theory  of  Markov  Chains,  the  steady  state  prob¬ 
abilities  of  finding  different  numbers  of  customers  at  the  facil¬ 
ity  at  departure  epochs  can  be  found  through  the  one-step  transi¬ 
tion  probabilities  of  this  Markov  chain.  It  involves  calculation 
of  the  transition  probabilities  and  solving  of  a  set  of  simultane¬ 
ous  equations,  the  total  number  of  which  depends  upon  the  total 
number  of  the  states  of  the  Markov  chain.  The  time  average  prob¬ 
abilities  of  finding  different  numbers  of  customers  at  the  service 
facility  are  then  related  to  the  steady  state  departure  probabili¬ 
ties  in  those  models  in  which  they  are  not  equal. 

Calculation  of  transition  probabilities  which  are  the  elements 
of  the  matrix  corresponding  to  the  set  of  simultaneous  equations. 


10 


may  be  difficult  in  complex  models.  But  the  number  of  equations 
to  be  solved  is  much  smaller  than  the  case  of  the  global  balance 
equations  in  the  generalized  method  of  stages  because  of  the 
smaller  state  space  of  the  Markov  chain. 

1.3.3  Supplementary  Variable  Technique 

A  method  of  making  a  non-Markovian  process  Markovian  is  to 
incorporate  the  missing  information  by  adding  a  continuous  vari¬ 
able  as  a  supplementary  variable  in  the  state  description.  In 
case  of  models  with  exponential  inter-arrival  times  and  arbitrary 
service  times,  the  state  of  the  system  is  defined  by  the  pair 
consisting  of  the  number  of  customers  at  the  facility  and  the 
elapsed  service  time  of  the  customer  already  under  service.  This 
approach  was  first  suggested  by  Kendall  [KEND  53],  but  was  first 
used  by  Cox  [COX  55a].  The  remaining  service  time  can  also  be 
used  as  the  supplementary  variable,  insteady  of  the  elapsed  ser¬ 
vice  time  [HEND  72].  As  compared  to  this  method,  a  discrete  vari¬ 
able,  namely  the  stage  of  the  service  time  distribution  of  the 
customer  already  under  service,  is  added  as  a  supplementary  vari¬ 
able  in  the  method  of  generalized  stages  described  in  1.3.1. 

After  the  supplementary  variables  are  included  in  the  state 
description,  balance  equations  can  be  written  and  solved.  The 
solution  of  the  equations  involves  transforms,  which  can  become 
very  difficult  for  complex  models. 
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1.4  REVIEW  OF  RELATED  WORK 


Analysis  of  single  class  finite  source  models  with  exponen¬ 
tial  inter-arrival  and  service  times  using  global  balance  equa¬ 
tions  can  be  found  in  books  on  queueing  theory  [GROS  74  and 
KLEI  75].  The  single  class  finite  source  model  with  exponential 
service  times  and  arbitrary  service  time  distributions  was  studied 
by  Jaiswal  {JAIS  68]  using  the  supplementary  variable  technique. 

In  this  model  the  proportion  of  departures  that  leave  a  certain 
number  of  customers  at  the  service  facility  is  not  equal  to  the 
proportion  of  time  the  same  number  of  customers  are  at  the  facil¬ 
ity.  Courtois  and  Georges  [COUR  71]  obtained  relations  between 
these  departure  average  and  time  average  probabilities  for  a 
M/G/l  queueing  model  with  state-dependent  arrival  and  service  pro¬ 
cesses  and  a  single  class  of  customers.  The  single  class  finite 
source  model  is  a  special  case  of  the  model  analyzed  by  them. 

Thir uvengadam  [THIR  65]  analyzed  a  M/G/l/  /N^.N^/PR  model 
with  non- preemptive  priority  rule,  using  the  supplementary  vari¬ 
able  technique  and  solved  the  resulting  differential  difference 
equations.  Jaiswal  and  Thir uvengadam  [JAIS  67]  considered  this 
method  tedious  and  modified  the  analysis,  basing  it  upon  the 
basic  server  sojourn  process,  which  starts  when  a  lower  priority 
customer  enters  service  and  ends  when  the  server  becomes  free 
to  accept  the  next  lower  priority  customer  for  service.  Using 
this  approach,  they  studied  a  M/G/l/  /N^jNj/PR  model  with  pre- 


emptive  priority  rule.  For  the  non-preemptive  service  rule,  they 
suggested  modifications  in  this  approach.  Extending  this  work, 
Jaiswal  [JAIS  68]  did  extensive  work  in  different  types  of  pri¬ 
ority  queueing  models.  For  a  M/G/l/  /N^j^/PR  model  with  non- 
preemptive  priority  rule,  he  derived  the  expected  length  of  the 
busy  period  duration,  in  terms  of  lengthy  functions  of  Laplace 
Transforms  of  service  time  distributions  of  the  classes,  using 
the  supplementary  variable  technique.  He  also  obtained  Laplace 
Transforms  of  the  joint  queue  length  probabilities  and  the  occu¬ 
pation  time  density.  These  are  difficult  to  solve  even  for  two 
classes.  Generalization  to  more  than  two  classes  presents 
immense  algebraic  difficulties.  Therefore,  Jaiswal  did  not  ex¬ 
tend  his  approach  to  more  than  two  classes.  We-Min  Chow  [CHOW  75] 
investigated  the  behavior  of  the  same  model  using  imbedded  Markov 
chain  analysis.  He  derived  a  relation  between  time  average  prob¬ 
abilities  and  departure  average  probabilities  in  terms  of  the 
integrals  of  functions  of  inter-arrival  and  service  time  distri¬ 
butions  of  the  customers  and  the  one  step  transition  probabilities. 
These  are  difficult  to  evaluate  even  for  two  classes. 

Benson  and  Cox  [BENS  51]  studied  a  single  finite  source 
priority  model  in  which  a  group  of  machines  failed  because  of  two 
types  of  failures.  Assuming  exponential  running  and  service  times, 
they  used  balance  equations  and  obtained  the  machine  availability 
and  the  operator  efficiency.  Jaiswal  and  Thiruvengadam  [JAIS  63] 
analyzed  a  similar  model  with  arbitrary  service  time  distribution 
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using  the  supplementary  variable  technique.  They  obtained  the 
steady  state  probabilities  of  the  number  of  machines  waiting  for 
repair,  the  operator  efficiency  and  the  machine  availability. 

None  of  the  above  works,  which  studied  either  the  multiple 
or  the  single  finite  source  priority  models  with  exponential 
inter-arrival  times  and  arbitrary  service  times,  could  obtain 
computationally  feasible  solutions  for  finding  the  system  per¬ 
formance  measures.  Most  of  these  studies  described  the  approach¬ 
es  which  could  be  used  to  analyze  such  models. 

Any  finite  source  queueing  model  can  be  represented  by  an 
equivalent  closed  network  with  two  nodes.  Queueing  networks 
play  an  important  role  as  performance  models  of  computers  and 
communications  systems.  Therefore,  considerable  research  has 
been  done  in  recent  years  on  queueing  networks .  This  has  re¬ 
sulted  in  a  significant  number  of  useful  results.  An  important 
result  among  these  is  the  discovery  of  Jackson  [JACK  63],  Gordon 
and  Newell  [GORD  67],  and  Basket,  et  al.  [BASK  75]  that  for  cer¬ 
tain  classes  of  networks  the  solutions  of  the  balance  equations 
are  in  the  form  of  a  product  of  simple  terms.  These  are  known 
as  product  form  solutions.  The  advantage  of  such  form  of 
solutions  is  that  the  problem  of  obtaining  the  probabilities 
reduces  to  that  of  normalizing  the  product  terms  to  form  a  pro¬ 
per  probability  distribution  and  of  the  computation  of  the 
normalizing  constant.  The  process  of  these  solutions  was  helped 


by  the  algorithms  published  by  Buzen  [BUZE  73]  and  Reiser  and 
Kobayashi  [REIS  75a,  REIS  75b,  and  REIS  77].  In  some  cases  it 
may  be  necessary  to  find  only  the  mean  values  of  the  queue  sizes, 
waiting  times,  and  utilizations  and  not  the  probability  distribu¬ 
tions.  If  in  such  cases  the  models  have  product  form  solution, 
then  a  solution  method  called  mean  value  analysis,  developed  by 
Reiser  and  Lavenberg  [REIS  80],can  be  used.  This  method  is  based 
on  the  relation  between  the  mean  waiting  time  and  the  mean  queue 
size  of  a  model  with  one  less  customer.  Without  the  necessity 
of  computing  the  normalization  constants,  this  technique  solves 
the  required  set  of  equations  numerically.  This  analysis  is 
considered  to  be  simpler  and  numerically  less  troublesome. 

There  are  many  models,  including  the  multiple  class  closed 
network  model  with  arbitrary  service  times  and  fixed  priority 
disciplines,  which  do  not  have  product  form  solution.  One  of 
the  ways  to  analyze  such  models  is  the  standard  method  of  writing 
the  global  balance  equations  and  solving  them,  if  service  times 
have  rational  Laplace  transforms.  For  large  models,  approximate 
numerical  solution  is  a  feasible  alternative.  Sauer  and  Chandy 
[  SAUE  80]  give  an  account  of  the  approximate  techniques  developed 
to  analyze  queueing  networks. 

The  approximate  techniques  are  basically  of  three  types. 

The  first  one  is  aggregation  which  replaces  subnetworks  by  com¬ 
posite  queues  which  will  produce  approximately  the  same  flow  of 


customers  through  the  queue  as  the  subnetwork  for  all  classes 


of  customers.  This  is  done  until  the  resulting  network  has  a 
feasible  solution.  This  was  applied  to  analyze  a  multiple  class 
closed  queueing  network  model  with  fixed  priority  discipline  by 
Sauer  and  Chandy  [SAUE  75].  The  second  type  is  diffusion  ap¬ 
proximation.  In  the  diffusion  approximation  the  variances  of 
the  inter-arrival  times  and  the  service  times  can  be  incorporated 
and  the  discrete  process  is  represented  by  a  continuous-time 
continuous-state  Markov  process.  Since  diffusion  processes  in 
complex  models  are  difficult  to  solve  algebraically,  heuristics 
are  used  to  simplify  the  solution,  thus  obtaining  a  diffusion 
approximation.  Diffusion  approximations  have  been  primarily 
used  to  open  networks  with  homogeneous  customers.  The  works  by 
Gaver  [GAVE  68]  and  Kobayashi  [KOBA  74]  are  some  examples  of 
this  approximation.  But  as  Kobayashi  [KOBA  78]  points  out, 
there  is  no  general  formula  available  that  helps  one  to  assess 
the  accuracies  of  the  solutions  obtained  through  the  approxima¬ 
tions  using  either  an  aggregation  or  a  diffusion. 

The  third  method  of  approximation  involves  using  mean 
value  analysis  in  the  case  of  models  which  do  not  have  product 
form  solutions.  Notable  work  in  this  direction  was  done  by 
Reiser  and  Lavenberg  [REIS  80],  Bard  [BARD  78],  and  Schweitzer 
[SCHW  79]  for  the  analysis  of  multiclass  closed  queueing  net¬ 
works.  But  the  results  are  generally  correct  only  in  the 
asymptotic  case.  Though  Bard  and  Schweitzer  claim  that  this 
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method  of  approximation  can  be  used  in  the  case  of  multiple 
class  networks  with  fixed  priority  discipline  and  non-exponen¬ 
tial  service  time  distributions,  no  such  known  analysis  is 
available. 

It  can  be  observed  from  the  preceding  review  that  no  feasi¬ 
ble  analytical  or  approximate  results  are  available  for  the 
multiple  finite  source  priority  model  with  non-preemptive  fixed 
priority  discipline  and  arbitrary  service  time  distributions. 
Simulation  is  the  most  widely  used  technique  to  study  this 
model.  But  simulation  is  expensive  and  time  consuming,  especial¬ 
ly  when  the  run  lengths  are  long  to  improve  the  accuracy  of  the 
results.  As  an  alternative,  recently  there  has  been  a  growing 
interest  in  numerical  algorithmic  methods  to  solve  the  models 
for  which  no  tractable  analytic  results  are  available.  One  of 
the  main  approaches  used  in  algorithmic  methods  is  to  represent 
the  state  probabilities  of  a  Markov  process  in  the  form  of  a 
set  of  linear  equations  and  to  find  an  efficient  solution  to 
these  equations  by  exploiting  the  special  structure  of  the  matrix 
corresponding  to  the  equations .  One  such  related  work  is  by 
Raju  and  Bhat  [RAJU  77]  who  developed  numerical  algorithms  to 
analyze  finite  waiting  room  capacity  queueing  models,  with  mul¬ 
tiple  classes. 
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1.5  OVERVIEW  OF  THIS  RESEARCH 


Chapter  2  contains  the  details  of  the  solution  procedure 
for  obtaining  the  mean  values  of  system  performance  measures 
which  lead  to  a  set  of  algorithms.  Based  on  the  logic  of  its 
development,  the  procedure  is  divided  into  parts  and  described 
in  separate  sections.  The  mathematical  relations  required  for 
computation  are  derived.  Then  the  algorithms  are  summarized  in 
a  way  suitable  for  adaption  as  a  computer  code  in  section  2.6. 
The  results  of  verification  of  the  algorithms  using  simulation 
are  illustrated  at  the  end  of  the  chapter  along  with  comments 
on  the  computational  aspects  of  the  algorithms. 

In  Chapter  3,  the  algorithms  developed  in  Chapter  2  are 
first  extended  to  obtain  time  average  marginal  and  joint  prob¬ 
abilities  of  the  number  of  customers  at  the  service  facility, 
using  Level  Crossing  Analysis.  Then  the  possibility  of  using 
the  algorithms  for  mixed  class  models  in  which  some  classes 
have  infinite  capacity  and  the  other  classes  have  finite  cap¬ 
acity  sources  is  discussed.  Finally,  conclusions  and  sugges¬ 
tions  for  further  research  are  given  at  the  end  of  the  chapter. 


& 
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CHAPTER  2 


SOLUTION  PROCEDURE 


2.1  INTRODUCTION 

Imbedded  Markov  chain  is  used  as  the  basic  method  of 
analysis  of  the  model  in  this  research.  Out  of  the  possible 
methods  of  approaches  discussed  in  Chapter  1,  the  supplemen¬ 
tary  variable  technique  is  found  to  be  a  very  difficult  method 
of  analysis  for  this  model,  especially  when  the  number  of 
classes  is  greater  than  2  [ JAIS  68] .  The  method  of  generalized 
stages  is  restricted  to  service  time  distributions  with  rational 
Laplace  transforms.  Even  in  cases  where  it  can  be  used,  the 
number  of  states  becomes  very  large  compared  to  that  of  im¬ 
bedded  Markov  chain. 

In  order  to  achieve  the  objective  of  this  research,  ef¬ 
forts  were  made  to  obtain  direct  formulae  for  the  mean  values 
of  system  performance  measures  which  can  be  related  to  the 
results  obtained  through  imbedded  Markov  chain  analysis.  This 
resulted  in  a  series  of  logical  and  sequential  steps  which  are 
described  in  the  following  sections. 
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2.2  BASIC  RELATIONS 


2.2.1  Relations  of  Mean  Values 

A  single  class  queueing  system  consisting  of  a  finite 
source  with  N  customers  served  by  a  single  server  can  be  re¬ 
presented  as  in  Figure  2.1.  Each  customer,  after  spending  an 
average  of  1/A  time  units  in  the  source,  arrives  at  the  service 
facility  and  waits  in  the  queue  for  time  units  on  the  average. 
The  mean  service  time  is  1/p  time  units  and  after  the  service  is 
completed,  the  customer  returns  back  to  the  source  without  any 
time  loss.  This  cycle  is  repeated  for  each  customer.  It  is 
assumed  that  the  system  is  in  steady  state  which  ensures  the 
existence  of  well  defined  and  finite  limiting  average  values  of 
1/A,  W  and  1/p  .  After  entering  the  system  marked  as  Box  A, 

q 

the  customer  never  leaves  it  until  departure  at  point  b  . 
Therefore,  Little's  formula  can  be  applied  to  the  whole  system 
[STID  78,  KOBA  78]  which  states  that  the  average  number  of  cus¬ 
tomers  within  the  system  (i.e.,  within  Box  A)  equals  the  average 
output  of  the  system  (at  point  b)  times  the  average  time  spent 
by  a  customer  within  the  system. 

The  average  number  of  customers  at  any  time  in  the  system 
is  equal  to  the  total  number  of  customers  in  the  system,  which 
is  equal  to  N,  as  the  customers  return  back  from  the  service 
facility  to  the  source,  without  any  time  loss.  The  average 
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output  of  the  system  is  equal  to  pp  where  p  is  the  utilization 
of  the  server.  This  is  valid  under  all  conditions.  The  only 
restriction  is  that  the  service  discipline  is  work  conserving 
in  that  the  service  times  of  the  individual  customers  are  not 
affected  by  the  discipline  [KOBA  78].  Therefore, 

N  =  pp(l/A  +  W  +  1/p) 

q 

which  gives 

Wq  =  N/pp  -  1/p  -  1/A  .  (2.1) 

No  assumption  was  made  about  the  distribution  of  either  the 
inter-arrival  times  or  the  service  times  in  deriving  relation 
(2.1).  The  only  requirements  are  that  the  system  should  be 
operating  under  steady  state  and  that  the  service  times  of  the 
individual  customers  are  not  affected  by  the  service  discipline. 

In  the  case  of  a  finite  population  model  with  capacity  N,  the 
effective  arrival  rate  of  the  customers  at  the  facility  is 
(N-L)A,  where  L  is  the  mean  number  of  customers  present  at  the 
facility  and  1/A  is  the  mean  inter-arrival  time  of  a  customer 
[GROS  74].  The  mean  amount  of  time  a  customer  spends  at  the 
facility,  W,  can  be  related  to  L  using  the  relation  L  =  (N-L)A  W. 

Using  this,  the  relation  W  =  +  1/p  and  (2.1),  L  can  be  ex¬ 

pressed  as, 
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•T 


L  =  N  -  pp/X  .  (2.2) 

L^,  the  mean  number  of  customers  waiting  in  the  queue,  is  re¬ 
lated  to  L  by  the  relation  L  =  L  +  p. 

q 

In  the  case  of  a  multiple  finite  source  system  which  has 
K(>^  2)  number  of  independent  sources  and  a  single  server,  the 
flow  of  customers  of  each  source  can  be  considered  separately 
in  a  subsystem  and  separate  boxes  like  Box  A  in  Figure  2.1  can 
be  constructed  for  each  subsystem.  If  the  service  discipline 
is  of  non-preemptive  fixed  priority  type,  which  is  work-conserv¬ 
ing  in  the  sense  described  earlier,  and  if  steady  state  condi¬ 
tions  can  be  assumed  to  exist,  then  Little's  formula  can  be  ap¬ 
plied  to  each  subsystem  as  before  and  the  following  relations 
are  obtained  for  i  *  1,2,...,K: 

(i)  Mean  waiting  time  of  class  i  customer  in  the  queue: 

Wqt  =  N./p.p.  -  l/v±  -  1/Xi  (2.3) 

(ii)  Mean  time  spent  by  a  class  i  customer  at  the  facility: 


W  =  W  +  1/p.  (2.4) 

1  qi  1 

(iii)  Mean  number  of  class  i  customers  present  at  the  facility: 


Li  =  Nt  -  Piui/Xi  ,  and  (2.5) 


(iv)  Mean  number  of  class  i  customers  waiting  in  the  queue: 


L  =  L.  -  p.  (2.6) 

qA  i  1 


where  refers  to  the  proportion  of  time  the  server  is  busy  with 
class  i  customers  or  the  time  average  probability  that  the  server 
is  busy  with  class  i  customers.  The  utilization  of  the  server  is 
given  by. 


K 

P  =  £P,  (2.7) 

i=l  1 

Relations  (2.3)  to  (2.7)  illustrate  that  the  mean  values  of 
the  system  performance  measures  of  interest  can  be  obtained  if 
the  p^'s  (i  =  1,2,...,K)  can  be  found.  In  the  next  subsection, 
the  method  of  obtaining  the  values  of  p^'s  is  described  using  the 
theory  of  regenerative  processes. 

2.2.2  Regenerative  Process 

A  regenerative  process  {X(t),  t  >  0}  is  a  stochastic  process 
which  starts  anew  probabilistically  at  an  increasing  sequence, 

0  <  Rx  <  R2  <  R3  ...,  of  random  epochs  on  the  time  axis  [0,®). 

Thus,  between  any  consecutive  epochs  and  R^+^,  the  portion 
{X(t),  <_  t  <  R^+^}  is  an  independent  and  identically  distri¬ 

buted  replicate  of  the  portion  between  any  other  two  consecutive 
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epochs  [ROSS  70].  Also,  the  time  interval  between  two  consecu¬ 
tive  epochs  R^  and  is  called  the  regenerative  cycle  i, 

whose  length  is  represented  by  T  .  These  lengths  of  the  regen¬ 
erative  cycles  are  independent  and  identically  distributed 
random  variables. 

Suppose  that  Y^  represents  a  reward  earned  during  the  re¬ 
generative  cycle  i  and  that  the  pairs  (T^ ,  Y  ) ,  SL  =  1,2 . are 

independent  and  identically  distributed.  If  Y(t)  denotes  the 
total  reward  earned  by  time  t,  then  the  limiting  value  of  the 
average  return  is  given  by  the  following  theorem. 

Theorem  2.1:  If  E(|Y^|)  and  ECT^)  are  finite,  then 

(i)  with  probability  1, 


Y(t)  E(Y) 
t  E(T) 


as  t  ->•  <*> 


(2.8) 


(ii) 


E(Y(t) )  ^  E(Y) 
t  E(T) 


t  -*  oo 


The  proof  of  this  theorem  can  be  found  in  [ROSS  70] .  According 
to  this  theorem,  the  expected  long-run  return  is  just  the  expected 
return  earned  during  a  cycle  divided  by  the  expected  length  of  a 
cycle. 

In  the  model  being  studied  in  this  research,  X(t)  can  re¬ 
present  the  total  number  of  customers  present  at  the  facility 
at  time  t  and  each  busy  cycle,  consisting  of  a  busy  period  and 
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an  idle  period,  can  be  considered  as  a  regenerative  cycle  [CHOW 


75].  Let  Y^  represent  the  amount  of  time  the  server  is  busy 

with  class  i  customers  (i  =  1,2,...,K)  during  the  busy  cycle  1 

and  Y(t)  denote  the  total  amount  of  time  the  server  is  busy 

with  class  i  customers  during  a  time  period  t.  Then  as  per 

Y(t) 

Theorem  2.1,  (which  is  equal  to  — —  as  t  -*•  ®)  can  be  ob¬ 
tained  by  taking  the  ratio  of  the  expected  length  of  time  that 
the  server  is  busy  with  class  i  customers  during  a  busy  cycle 
to  the  mean  length  of  one  busy  cycle.  Let  E(B)  be  the  expect¬ 
ed  length  of  one  busy  period,  ECB^) ,  (i  =  1,2,...,K),  be  the 
expected  length  of  the  busy  period  during  which  the  server  is 
busy  with  class  i  customers,  and  E ( I )  be  the  expected  length 
of  one  idle  period.  Then,  as  per  Theorem  2.1,  for  i  =  1,2,...,K, 


E(Bi) 

Pi  =  E(B)  +  E(I) 


(2.9) 


where 


K 

E(B)  =  L  E(B. ) . 
i=l  \ 


The  idle  period  starts  at  the  instant  when  all  the  cus¬ 
tomers  are  at  their  respective  sources.  Since  the  inter-arrival 
time  of  each  customer  of  class  i,  (i  *  1,2,...,K),  follows  an 
exponential  distribution  with  mean  1 as  per  the  model  des- 
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cription,  the  mean  time  interval  after  which  the  first  cus¬ 


tomer  of  class  i  arrives  at  the  service  facility,  after  the 
start  of  an  idle  period  is  1/N^X^.  A  customer  of  any  class 
can  terminate  an  idle  period  and  therefore  the  mean  length  of 
the  idle  period  is  given  by 

K  -1 

E(I)  =  [  Z  NX]  1  .  (2.10) 

i=l 

Since  in  relation  (2.9)  only  the  values  of  E(B^)’s  are 
unknown,  the  next  step  in  the  solution  development  is  to  find 
the  values  of  EtB^),  (i  =  1,2,...,K).  Imbedded  Markov  chain 
analysis  is  used  for  this  purpose. 

2.2.3  Imbedded  Markov  Chain  Analysis 

As  discussed  in  Chapter  1,  in  the  model  being  studied  in 
this  research,  the  epochs  at  which  customers  depart  from  the 
service  facility  constitute  a  set  of  regenerative  points.  As 
this  model  has  K  distinct  classes  of  customers,  the  numbers  of 
customers  of  these  K  different  classes  present  at  the  service 
facility  at  departure  epochs  form  a  Markov  process.  The  states 
of  the  Markov  chain  of  this  process  are  expressed  by  vectors 
which  consist  of  K  elements  corresponding  to  the  numbers  of 
customers  of  K  classes  present  at  the  facility  at  the  departure 
epochs.  Let  X  and  X  , .  denote  the  states  of  the  process  at 


the  nth  and  (n+l)St  departure  epochs,  respectively,  (n  =  1,2,...). 

Let  the  elements  of  the  row  vectors  =  (j^,  j  j  *  •  •  • » j  ^  •  t  jg) 

=  (u^,U2 , . . .  ,u^ , . . .  ,Uj,)  represent  the  numbers  of  customers  of 

ttl 

the  corresponding  classes  present  at  the  facility  at  the  n  and 
s  t 

the  (n+1)  departure  epochs,  respectively.  Then  the  elements 

of  the  transition  probability  matrix  corresponding  to  this  Markov 

chain  are  given  by  the  conditional  probabilities,  P[X^+^  =  1L |X  =  j], 

for  different  possible  values  of  the  elements  of  u^  and  j_.  The 

total  number  of  states  of  the  Markov  chain  is  given  by  the  number 

of  all  possible  combinations  of  the  numbers  of  customers  of  K 

classes  that  can  be  left  behind  at  the  facility  by  a  departing 

K 

customer.  This  is  equal  to  {  Z  (N.+l)}  -  1,  which  is  denoted  as 

i=l  1 

(m+1)  for  the  sake  of  simplicity. 

As  this  type  of  imbedded  i.arkov  chain  is  finite,  aperiodic, 
and  irreducible,  all  its  states  are  ergodic  [KLEI  75].  This  en¬ 
sures  the  existence  of  steady  state  probabilities.  Let  the  row 
vector  containing  the  steady  state  probabilities  of  this  Markov 
chain  be  represented  by  A.  It  has  (m+1)  elements  corresponding 
to  the  (m+1)  states  of  the  Markov  chain.  Two  ways  are  used  to 


represent  the  elements  of  this  vector.  One  way  is  to  represent 


a  typical  element  of  A  as  A^ ,  which  refers  to  the  j  element 
of  A.  In  the  second  way  of  representation,  A(v)  refers  to  the 


steady  state  probability  that  just  after  the  departure  of  a 
customer,  the  state  of  the  process  is  v  =  (v^ ,v^ , . . • ,v^ , . . . ,v^) . 
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The  second  way  of  representation  is  used  whenever  the  informa¬ 
tion  regarding  the  number  of  customers  left  behind  by  a  depart¬ 
ing  customer  is  of  importance. 

2. 2. 3.1  The  Relation  Between  E(B^)  and  A 

Let  E(R)  represent  the  expected  number  of  customers  served 
by  the  server  during  a  busy  cycle  and  p^  represent  the  propor¬ 
tion  of  class  i  customers  served  during  a  busy  cycle,  i.e.,  the 
steady  state  probability  that  a  departing  customer  is  of  class 
i.  Then  EfR^),  the  mean  number  of  class  i  customers  served  dur¬ 
ing  a  busy  cycle,  is  given  by,  for  i  *  1,2 . K, 

E(Ri)  =  E(R)  p±  (2.11) 


and  therefore 


E(B.)  =  E(R1)/lji 

=  E(R)  Pj/p^  .  (2.12) 

As  per  the  second  way  of  representing  the  elements  of  A, 

A(0)  stands  for  the  steady  state  probability  that  a  departing 
customer  leaves  behind  an  empty  facility,  where  0  stands  for 
the  vector  with  all  zero  elements,  i.e.,  £  *  (0,0,...,0).  This 
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implies  that  A(0)  is  the  steady  state  time  average  proportion 
of  served  customers  who  terminate  busy  periods,  and  that  on  the 
average,  one  out  of  every  [A(0)]  ^  customers  terminate  a  busy 
period.  Therefore, 


E(R)  =  [ A (0) ] 


(2.13) 


The  quantity  can  be  divided  into  two  parts:  p^  and  p^* 

The  first  part,  p  is  the  steady  state  probability  that  immedi¬ 

ately  after  the  termination  of  an  idle  period,  the  first  depart¬ 
ing  customer  is  of  class  i.  As  the  steady  state  probability 
that  a  departing  customer  starts  an  idle  period  is  A(0)  and  the 

probability  that  an  idle  period  is  terminated  by  a  class  i  cus- 
K  -1 

tomer  is  l  N^X^]  »  (JAIS  68],  p^  is  given  by 

i=l 

K  -1 

Pil  =  A(0)  NiXi(  I  NtAt]  A  (2.14) 


The  second  part,  P^.  is  the  steady  state  probability  thst 
excluding  the  first  departing  customer  during  a  busy  period,  any 
other  departing  customer  is  of  class  i.  If  the  elements  of  the 
row  vector  v  =  (v^,v2» • • • »v^, . . . ,v^)  represents  the  number  of  cus¬ 
tomers  of  the  corresponding  classes  left  behind  by  a  departing 
customer,  then  a  class  i  customer  will  be  the  next  customer  to 


be  served  and  to  depart  as  long  as  v^  )*  0  and  Vj  =  0  for  all 
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j  <  i.  This  is  because  of  the  fixed  priority  service  discipline. 
Therefore, 


p.,  =  2  [  I  {  r  A ( v) } ] 

v.=l  £=i+l  v  =0 

1  a 


(2.15) 


where  =  0  for  all  j  <  i,  in  As  p^ 
(2.11)  to  (2.14),  for  i  =  1,2,...,K, 


p^  +  and  from 


K 


-1 


E(B.)  =  l/p^N  A  (  Z  HJ)  *  + 
i=l 


N. 

l 


N, 


K  l 
A(0)  Z  (  Z  {  Z  A(v)  })  1 
v_L=l  lri+1  v  £=0 


(2.16) 


where  vj  =  0  f°r  all  j  <  *•»  in  v. 

In  equation  (2.16),  the  only  unknown  quantities  are  the 
elements  of  the  steady  state  probability  vector,  A.  Therefore, 
the  next  step  in  the  solution  procedure  is  to  find  the  values  of 
the  elements  of  A. 


2. 2. 3.2  Obtaining  The  Elements  of  A 

Let  P  represent  the  transition  probability  matrix  of  the 
Markov  chain.  Then  the  vector  A  can  be  uniquely  determined  by 
solving  the  system  of  linear  homogeneous  equations 

A(I-P)  =  0  (2.17) 
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subject  to  the  normalizing  condition 


nri-L 

I  A  =  1  •  (2.18) 

j=l  J 

In  (2.17),  I  represents  the  identity  matrix  of  size  (nrfl)x(nH-l) . 

Of  the  (nrt-1)  equations  in  (2.17),  only  m  are  independent.  So 

the  matrix  (I-P)  does  not  have  an  inverse.  If  one  of  the  (nri-1) 

equations  in  (2.17)  is  eliminated  then  the  remaining  m  equations 

along  with  (2.18)  can  be  used  to  obtain  a  unique  solution  for  A. 

Let  the  first  equation  in  (2.17)  be  removed,  which  is  equivalent 

to  eliminating  the  first  column  of  (I-P).  Expressing 

A  in  terms  of  A  . , ,  the  system  of  equations  in  (2.17)  can  be 
m  nrrjL 

written  as 


aVf)1  -  Vl  *  • 


(2.19) 


In  (2.19)  the.  row  vector  A^  contains  the  first  m  elements 
of  A,  the  matrix  (I-P)^  is  of  size  mxm  and  consists  of  all  elements 
of  (I-P)  except  those  in  the  last  row  and  the  first  column  of  (I-P) 
and  the  row  vector  Z  contains  the  negative  of  the  last  m  elements 
of  the  last  row  of  (I-P).  Now  the  m  equations  in  (2.19)  are  in¬ 
dependent  and  therefore  the  matrix  (I-P)^  has  an  inverse.  Multi¬ 
plying  both  sides  of  (2.19)  by  the  inverse  of  (I-P)\  the  elements 
of  A^  are  given  by. 


! 


A1  =  A^UHI-P)1]'1}  . 

By  using  equations  (2.20)  and  (2.18),  the  values  of  all  the 
steady  state  probabilities  can  be  obtained. 

At  this  stage,  the  problem  reduces  to  the  following: 

(i)  Generation  of  the  elements  of  the  transition  prob¬ 
ability  matrix  P  and  from  these,  obtaining  the 
elements  of  the  row  vector  Z  and  of  the  matrix 
(i-P)1; 

(ii)  Inversion  of  (I-P)^"; 

(iii)  Obtaining  the  values  of  A^ ,  j  =  l,2,...,m+l, 
using  (2.20)  and  (2.18). 

These  three  steps  are  described  in  the  following  sections. 


2.3  GENERATION  OF  THE  ELEMENTS  OF  P 

As  stated  in  section  2.2,  the  elements  of  P  are  given  by 
the  conditional  transition  probabilities,  **  ujx^  -  Jl]  . 

The  calculation  of  these  probabilities  is  discussed  under  two 
sets  of  conditions,  depending  upon  whether  the  n^  departing 
customer  leaves  behind  a  non-empty  or  an  empty  facility. 

Under  the  first  set  of  conditions,  the  elements  of  J_,  which 
gives  the  state  of  the  process  at  the  nC^  departure  epoch,  are 
such  that  there  is  at  least  one  non-zero  element,  j^,  (£  ■  1,2, 
...,K),  with  j  =  0  for  all  i  <  l,  if  l  4  1.  This  implies  that 
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the  n  departing  customer  leaves  behind  at  least  one  customer 
of  class  l  at  the  facility  without  any  customers  of  higher  pri¬ 
ority  classes  waiting  for  service.  Then,  because  of  the  fixed 

3 1 

priority  discipline,  the  next,  i.e.,  the  (n+1)  ,  customer  to  be 

serviced  belongs  to  class  l.  Therefore,  in  order  to  have  the 

s  t 

state  of  the  process  at  the  (n+1)  departure  epoch  to  be 

u  *  (u^.u^, • . . .u^) ,  the  number  of  arrivals  during  the  service 
s  t 

period  of  the  (n+1)  customer  has  to  be  (u^  -  j^+1)  for  class  J l 
and  (u^  -  j ^)  for  classes  i  =  1,2,...,K,  but  not  l.  Therefore, 
the  conditional  probability  is  given  by 

P[XR+1  *  uJXjj  “  jj  =  Probability  [number  of  arrivals  during 

the  (n+1) st  service  time 

■  *  (ui~ji) » • '  * » 

*Ujt-l'J  Jl-1^  *  ^UH+l“-U+l^  ’ 

• • • »  ^UK"JK)  1 3^  3  1  and  the 

(n+1) st  service  time  is  that 
of  a  class  l  customer,  i.e., 

V* 

As  per  the  description  of  the  model,  the  sources  are  finite 
and  independent  of  each  other  and  the  times  spent  by  the  customers 
of  each  class  at  the  corresponding  sources  are  exponentially  dis¬ 
tributed  random  variables.  Therefore,  the  number  of  arrivals  of 
each  class  customers  at  the  facility  within  a  given  time  inter¬ 
val  is  a  binomially  distributed  random  variable,  independent  of 
the  arrivals  of  other  class  customers.  Let  pOn^n^t)  denote  the 


probability  that  the  number  of  arrivals  of  class  i  customers, 

(i  =  1,2,...,K),  at  the  facility  within  a  time  period  t  is  m^ 
when  there  are  n^  customers  of  class  i  at  the  facility  at  the 
beginning  of  the  time  period.  Then, 

(N,-n,\  -X.t  tn  -X  t  N.-n.-m. 

1  1  (1-e  1  )  (e  1  )  1  1  1  (2.21) 

mi  / 

where  m^.n^  -  0,1,2 . N^ ,  mi+ni  i.  N^,  and  1/X^  is  the  mean 

time  spent  by  a  class  i  customer  at  the  source. 

Now,  for  ^  0  and  j.  =  0  if  i  <  i  and  l  f  1,  the  condi¬ 
tional  probability  can  be  written  as 

K 

p(x +1-«lxn"ii  =  A  n  p(ui"WsP  Jp(urJt+1;Ja»s£)dF(8P»  (2-22) 

o  i=l 

tt 

where  F^^)  is  the  distribution  function  of  the  service  time  of 
class  l  customer.  After  substituting  the  values  of  p(u^-j^;j^,s  j) 
and  P(u£-j£+l'» j^*^)  from  (2.21),  equation  (2.22)  becomes,  for 
1*  0  and  =  0  if  i  <  l. 


6i  o'  \  dF<si> 


where 


K 

{  n 

i=i 


Mi' 

iui_;ji 


V*i 

.vV1) 


(2.23) 


(2.24) 
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and 


K  ~Aisi  ui"Ji  “Aisi  Ni~ui  ~Xlal  ul~h+1  ■Xi.s£ 

0,  -{  n  (l-e  1  l)  1  i(e  1  l)  1  1}(l-e  1  1  *  (e  *  £)  *  1 

i-1 

H 


(2.25) 


Under  the  second  set  of  conditions,  the  nC^  departure  leaves 

behind  an  empty  facility,  i.e.,  in  vector  j  ■  0  for  i  -  1,2, 

s  t 

. . . ,K,  thereby  starting  an  idle  period.  Then  the  (n+1)  service 
time  depends  on  the  class  of  customer  who  terminates  this  idle 
period.  This  situation  is  illustrated  by  Figure  2.2  in  which 
the  n~k  customer  leaves  at  A  starting  an  idle  period.  The  cus¬ 
tomer  who  arrives  at  B  terminates  the  idle  period  and  after  being 


FIRST  ARRIVAL 


4  -  (0,0, . . . ,0) 


MORE  ARRIVALS 


X 

-n+1 


(uj ,U2 , • .  • , 


IDLE  PERIOD 


(n+l)St  SERVICE  TIME 


nth  DEPARTURE 


(n+l)8t  DEPARTURE 


BEHAVIOR  OF  THE  PROCESS  WHEN  j  -  0,  FOR  i  -  1,2,...,K. 

FIGURE  2.2 
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s  t 

served  becomes  the  (n+1)  departure  at  C.  From  B  to  C  the  pro¬ 
cess  behaves  as  if  at  point  B  there  is  a  departure  leaving  be¬ 
hind  only  that  customer  who  terminated  the  idle  period.  Because 
any  class  customer  can  terminate  the  idle  period  at  B,  the  re¬ 
quired  conditional  probability  is  given  by 


PI^i 


-  u  X 
— n 


-  0]  = 


K 

£  {Prob[idle  period  is  terminated  by  a  class 
i=l 

i  customer]  Prob[X  ,.=u|x  =1]} 

— n+j. - n 


(2.26) 


where  ()  refers  to  the  row  vector  containing  all  zero  elements  and 

the  elements  of  the  vector  ^  are  given  by  *  1  if  £  =  i  and 

j^»0ifi/i.  As  stated  in  section  2. 2. 3.1,  the  probability 

that  an  idle  period  is  terminated  by  a  class  i  customer  is  equal 
K  -1 

to  N.X  [  I  N  .1.]  .  So  equation  (2.26)  becomes 

1  1  il=l  11 


(2.27) 


where  P[X  =u|X  >*jj  with  j  **  1  if  l  =  i  and  j  ■  0  if  l  ^  i, 
n+1  ~n 

can  be  obtained  from  (2.23). 

The  evaluation  of  the  integral  in  (2.23)  depends  upon  F(s^) 

8 1 

which  is  the  distribution  function  of  the  (n+1)  service  time. 

In  Appendix  B  the  integral  is  evaluated  when  the  density  function 
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cf  is  exponential,  hypo-exponential,  or  hyper-exponential  for 
the  purpose  of  illustration-  These  are  the  density  functions 
frequently  used  to  represent  service  demands  [KOBA  78]  because 
these  density  functions  have  rational  Laplace  transforms.  In 
computer  systems,  the  coefficients  of  variation  (the  ratio  of  the 
standard  deviation  to  the  mean)  of  the  service  time  distributions 
are  found  to  be  greater  than  1  {ANDE  72]  and  thus  hyper-exponen¬ 
tial  density  function  is  appropriate  in  such  cases. 

Once  the  elements  of  the  transition  probability  matrix  are 
generated,  the  elements  of  (I-P)^  and  the  row  vector  Z  can  be 
easily  obtained.  The  next  step  then  is  to  find  the  inverse  of 
(I-P)^  which  is  described  in  the  next  section. 

2.4  INVERSION  OF  THE  MATRIX  (I-P)1 

2.4.1  Introduction 

The  inversion  of  a  matrix  becomes  tedious  and  time  consum¬ 
ing  as  its  size  increases.  The  size  of  (I-P)^  is  mxm  where  m  is 
K 

equal  to  {  II  (N,+l)}-2.  The  value  of  m  becomes  very  large  even 
i=l  1 

for  small  values  of  N^s.  One  of  the  goals  in  the  development 
of  the  solution  procedure  is,  therefore,  to  develop  an  efficient 
numerical  algorithm  for  inverting  (I-P)^. 

A  good  numerical  technique  used  for  inverting  any  matrix 
should  take  into  account  and  exploit  the  structure  of  that  matrix. 


i 
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In  some  cases,  it  may  be  possible  to  alter  the  structure  of  a 
matrix  so  as  to  suit  it  to  an  efficient  inversion  technique.  In 
the  case  of  the  matrix  (I-P)\  in  addition  to  altering  the  struc¬ 
ture  of  the  matrix,  it  is  necessary  also  to  modify  an  available 
inversion  technique  to  suit  the  altered  structure. 

There  are  some  interesting  matrix  structures  which  are  en¬ 
countered  in  the  case  of  finite  Markov  chains  of  queueing  models. 
Among  those,  the  two  types  most  related  to  this  research  are  al¬ 
most  left  triangular  and  left  triangular  structures  [RAJU  77]. 
Definition  1;  A  nxn  square  matrix  C  is  almost  left  triangular 
if  its  elements  are  such  that,  for  i  =  l,2,...,n, 

c.  .  =  0  for  all  j  >  i+1  . 

1 » J 

Definition  2:  A  nxn  square  matrix  C  is  left  triangular  if  its 
elements  are  such  that,  for  i  =  l,2,...,n, 

c,  .  =  0  for  all  j  >  i. 
i»J 

At  this  stage  it  is  necessary  to  introduce  the  terminologies,  the 
superdiagonal  and  the  main  diagonal  of  a  square  matrix. 

Definition  3:  the  super  diagonal  of  a  nxn  square  matrix  C  con¬ 
sists  of  the  elements  c,  ., 

*■  *  J 


such  that  j  =  i+1  for  i  *  l,2,...,n. 


Definition  4:  The  main  diagonal  of  a  nxn  square  matrix  C  con¬ 
sists  of  the  elements  c  ,  such  that  j  =  i  for  i  =  l,2,...,n. 

i  >  J 

If  the  structure  of  the  matrix  (I-P)  in  equation  (2.16)  is 
almost  left  triangular,  then  the  matrix  (I-P)1  has  a  left  tri¬ 
angular  structure.  The  advantage  of  having  a  matrix  with  left 
triangular  structure  and  all  non-zero  elements  along  the  main 
diagonal  is  that  computationally  efficient  methods  are  available 
for  the  inversion  of  such  a  matrix  (SCHR  73].  But,  in  the  case 
of  the  model  in  this  research,  it  is  not  possible  to  obtain  a 
left  triangular  structure  with  all  non-zero  main  diagonal  ele¬ 
ments  for  the  matrix  (I-P)1.  It  is,  however,  possible  to  place 
most  of  its  elements  into  a  left  triangular  form  while  having 
only  a  few  elements  above  the  main  diagonal  by  a  suitable 
arrangement  of  Its  column  and  row  states.  This  arrangement  is 
also  advantageous  because  the  inverse  of  the  left  triangular 
part  .of  (I-P)1  can  be  obtained  first  using  the  already  avail¬ 
able  technique  and  then  this  inverse  can  be  modified  taking  the 
elements  above  the  main  diagonal  into  consideration.  This  di¬ 
vides  the  procedure  of  inverting  (I-P)1  into  three  basic  parts. 
They  are,  (i)  arrangement  of  the  row  and  column  states  of  (I-P)1 
to  place  most  of  its  elements  into  a  left  triangular  form  leaving 
the  remaining  elements  above  the  main  diagonal;  (ii)  inversion 
of  the  lower  triangular  part  using  an  available  numerical  tech¬ 
nique;  (iii)  modification  of  the  inverse  obtained  in  (ii),  taking 
the  elements  above  the  main  diagonal  into  consideration. 


These  are  explained  in  detail  in  the  following  subsections. 


2.4.2  Arrangement  of  the  States 

The  arrangement  of  the  states  of  (I-P)^  is  done  prpgressive- 
ly,  starting  with  the  state's  of  the  transition  probability  matrix 
P,  then  the  states  of  the  matrix  (I-P),  and  finally  ending  up 
with  matrix  (I-P)^.  As  matrix  (I-P)'*'  is  obtained  by  deleting 
the  first  column  and  the  last  row  of  (I-P),  the  main  purpose  of 
the  arrangement  of  the  states  is  to  group  most  of  the  elements  of 
(I-P)  into  an  almost  left  triangular  form  with  non-zero  elements 
along  the  super  diagonal.  This  is  helped  by  the  following  facts: 

(i)  the  imbedded  Markov  chain  of  the  model  has  unit  jumps  at 
regeneration  points  which  means  that,  at  departure  epochs,  the 
number  of  customers  at  the  facility  decreases  by  at  most  one; 

(ii)  because  of  the  fixed  priority  rule,  the  number  of  customers 
of  class  l  (i  =  2,3,. ..,K)  can  decrease  by  one  at  departure  epochs 
only  if  there  are  no  customers  of  all  higher  priority  classes 
present  at  the  facility  at  the  previous  departure  epoch. 

2.4. 2.1  Arrangement  of  the  States  of  P 

A  typical  row  state  in  the  matrix  P  represents  the  number 
of  customers  left  behind  at  the  facility  by  the  n^  departing 
customer  (n  =  1,2,...).  It  is  denoted  as  *  j_,  where  j_  is  the 
row  vector  containing  the  elements  j^, j2» ... ,jK  which  represent 
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the  numbers  of  customers  of  corresponding  classes.  Let  X^ 

stand  for  the  row  state,  i  =  1,2 , . . . , (nri-1) .  A  typical 

column  state  in  P  represents  the  number  of  customers  left  be- 

s  t 

hind  by  the  next,  i.e.,  (n+1)  departing  customer.  It  is 

denoted  as  =  u  where  u  is  the  row  vector  containing  the 

elements  u.  , u_ , . . . ,u„  which  represent  the  numbers  of  customers 
1  /  K. 

v  th 

of  corresponding  classes.  Let  X^+^  stand  for  the  v  column 
state.  In  P  the  arrangement  of  row  states  is  identical  to 
that  of  the  column  states,  i.e.,  for  i  =  1,2 . m+1,  ^X  = 

V 

X  . , ,  when  v  =  i . 

—n+1 

The  method  chosen  to  arrange  the  states  in  P  forms  a  con¬ 
venient  basis  from  which  the  final  arrangement  of  the  states  of 
(I-P)  and  (I-P)^  emerge,  yielding  the  desired  matrix  structures. 
This  method  of  arrangement  is  explained  with  respect  to  the  row 
states  which  are  arranged  as  per  the  following  steps: 

(i)  Set  i  =  1  and  let  the  first  row  state  from  the 
top  be  equal  to 


x^  =  (n1,n2,...,nk_1,  Nk-1). 


Set  j£  =  N£  for  l  =  1,2,...,K-1  and  jR  =  NR-1. 

K 

(ii)  Set  i  =  i+1  and  U  =  K.  If  i  <  {  n  (N  +1)}-1, 

k=l 

go  to  step  (iii).  Otherwise,  go  to  step  (vi) . 
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(iii)  If  in  ft  0,  set  Ju 


j^-1  and  go  to  step  (v) . 


Otherwise,  go  to  step  (iv) . 

(iv)  Set  and  ill  •  8.1-1.  Go  to  step  (iii). 

(v)  Set  the  iC^  row  state 

2^  "  (Ji» ^2*  *  ’  ’  ' 

Go  to  step  (ii). 

(vi)  Stop. 


Because  the  maximum  value  of  i  is  set  as  m+1 
the  last  row  state  generated  is 


K 

{  n  (N,+l)  }-l, 
k-1  K 


^fl  -  (0,0 . 0). 

To  illustrate  the  above  steps,  an  example  is  considered  in 

which  there  are  three  classes  of  customers  with  ■  2,  N2  ■  1, 

and  ■  2.  As  per  step  (i)  the  first  state  is  ■  (2,1,1). 

2 

When  i  *=  2,  steps  (iii)  and  (v)  result  in  ■  (2,1,0).  When 
i  **  3,  as  ■  0,  step  (iv)  changes  to  ■  2  and  step  (iii) 

3 

changes  to  1-1  ■  0,  resulting  in  ■  (2,0,2).  The  next  two 
4  5 

states  are  X  -  (2,0,1)  and  X  -  (2,0,0).  When  i  ■  6,  because 
— n  — n 

and  j2  are  zeroes,  steps  (iii)  and  (iv)  change  to  ■  2, 
J2  to  ^2  "*  anc*  ^1  t0  Therefore  5C^  ■  (1,1,2).  Pro- 
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3 

ceeding  in  this  way,  when  i  =  {  II  (N, +1)  }-l  =  17,  the  last  row 

k=l  * 

17  K  x 

state  obtained  is  X  =  (0,0,0).  The  column  states  are  ordered 
— n 

in  the  same  way  as  the  row  states  and  the  resulting  P  matrix  is 
given  in  Figure  2.3. 

To  help  in  the  modification  of  the  arrangement  of  the  states 
in  (I-P)  matrix,  the  row  states  in  P  are  now  grouped  into  K  ordered 
row  sets  based  on  their  elements.  The  row  states  containing  the 
maximum  number  of  customers  of  the  first  (K-l)  classes  are  grouped 
together  to  form  row  set  1.  These  are  the  first  N  row  states 

iX 

from  the  top.  The  row  states  containing  the  maximum  number  of 
customers  of  first  (K-2)  classes  only  are  then  grouped  together 
as  row  set  2.  These  are  the  next  NK_1<N},+1)  row  states.  In 
general,  the  row  set  kl,  (kl  =  1,2,...,K),  contains  those  row 
states  with  the  maximum  number  of  customers  of  the  first  (K-kl) 
classes  only  and  the  total  number  of  row  states  in  this  set  is 
given  by 


r  (kl)  =  N. 


K-kl+1 


K 

H  (N  +1)} 
£=K-kl+2 


The  last  row  set  K  does  not  contain  the  maximum  number  of  cus¬ 
tomers  of  any  class.  The  column  states  of  P  are  also  grouped  in 
the  same  way  as  the  row  states,  starting  from  the  left.  This 
arrangement  of  the  states  does  not  place  most  of  the  elements 
of  P  in  an  almost  left  triangular  structure. 


;  Nz  -  1;  N3  -  2)  CS^  -  ORDERED  COLUMN  SET  1 

CS£  -  ORDERED  COLUMN  SET  2 
CS’  "  ORDERED  COLUMN  SET  3 
RSl  -  ORDERED  ROW  SET  X 
RS2  -  ORDERED  ROW  SET  2 
RS3  -  ORDERED  ROW  SET  3 

P  MATRIX 


FIGURE  2.3 


In  the  example  considered,  the  rows  are  divided  into  3 
ordered  sets.  The  row  set  1  consists  of  those  rows  which  con¬ 
tain  the  maximum  number  of  customers  of  classes  1  and  2;  and 
there  are  r(l)  =  2  of  such  rows.  The  next  rows  containing  the 
maximum  number  of  customers  of  class  1  only  are  grouped  as  row 
set  2  and  the  number  of  rows  in  this  set  is  given  by  r(2)  = 
1{(2+1)}  =  3.  The  third  row  set  consists  of  r(3)  =  2{ (1+1) (2+1) } 
*  12  rows  which  do  not  contain  the  maximum  number  of  any  class 
customers.  Column  states  are  also  grouped  in  the  same  way  and 
this  arrangement  is  illustrated  in  Figure  2.3. 

The  row  states  and  the  column  states  of  their  respective 
first  (K-l)  sets  contain  the  maximum  number  of  class  1  cus¬ 
tomers.  The  main  diagonal  elements  in  P  along  these  rows  and 
columns  are  all  zeros  as  the  corresponding  row  and  column  states 
cannot  communicate  with  each  other.  The  super  diagonal  of  P 
contains  zero  elements  corresponding  to  the  noncommunicating 
row  and  column  states. 

2. 4. 2. 2  Arrangement  of  the  States  of  (I-P) 

Matrix  (I-P)  is  obtained  by  subtracting  P  from  the  identity 
matrix  I  of  size  (m+l)x(nri-l).  All  the  non-zero  elements  of  P 


i 


except  the  main  diagonal  elements  become  negative  in  (I-P) .  The 
zero  elements  along  the  main  diagonal  of  P  become  1  and  all  other 
zero  elements  remain  as  zeros  in  (I-P).  This  implies  that  the 


main  diagonal  elements  of  the  first  (K-l)  row  sets  of  (I-P) 
are  1  and  the  super  diagonal  of  (I-P)  contains  some  zero  ele¬ 
ments.  Also  a  major  part  of  the  elements  of  (I-P)  do  not  form 
an  almost  left  triangular  structure. 

To  group  most  of  the  elements  of  (I-P)  into  a  left  triangu¬ 
lar  structure  and  to  have  all  non-zero  elements  along  the  super 
diagonal,  the  column  states  of  (I-P)  are  reordered  keeping  the 
ordering  of  row  states  unchanged.  The  arrangement  of  the  column 
states  can  be  different  from  that  of  the  row  states  in  (I-P) ,  be¬ 
cause  the  order  of  column  states  represents  only  the  order  in 
which  the  linear  equations  in  (2.16)  are  arranged.  The  reorder¬ 
ing  of  the  column  states  is  done  by  placing  the  columns  of  the 
last  column  set,  K,  of  P  in  the  first  positions  from  the  left  of 
(I-P),  the  columns  of  the  column  set  (K-l)  in  the  next  positions, 
and  so  on  and  placing  finally  the  columns  of  the  set  1  of  P  in 
the  last  positions  of  (I-P).  The  order  of  columns  within  a  set 
is  not  altered.  This  procedure  is  formally  stated  in  the  follow¬ 
ing  steps: 

(i)  Renumber  the  column  sets  in  P  in  such  a  way  that 
set  kl,  (kl  =  1,2,...,K)  becomes  set  k2  where 
k2  -  K-kl+1 . 

(ii)  Rearrange  the  sets  in  the  ascending  order  of  k2 
starting  from  the  left  in  (I-P). 
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Applying  these  steps  to  the  example,  column  set  3  of  P 
which  consists  of  the  12  column  states  from  (1,1,2)  to  (0,0,0) 
now  becomes  column  set  1  in  (I-P) ;  column  set  2  of  P  which  con¬ 
sists  of  3  column  states  from  (2,0,2)  to  (2,0,0)  now  becomes 
the  column  set  2  in  (1-P)  and  finally  the  column  set  1  becomes 
the  last  column  set  3  in  (I-P).  It  is  illustrated  in  Figure  2.4. 

This  reordering  of  the  column  states  places  all  the  elements 
of  (I-P),  except  the  unity  elements  along  the  rows  of  the  first 
(K-l)  row  sets,  in  an  almost  left  triangular  structure.  Because 
of  the  way  in  which  the  column  states  are  now  ordered,  these 
unity  elements  are  located  above  the  super  diagonal  in  column 
sets  2  through  K  and  arranged  along  (K-l)  lines  parallel  to  the 
super  diagonal.  If  these  lines  are  numbered  from  left  to  right, 
then  the  unity  elements  along  line  t  lie  in  column  set  (t+1)  and 
row  set  (K-t).  There  are  r(K-t)  unity  elements  along  line  t, 
where  r(K-t)  can  be  obtained  from  (2.28).  The  total  number  of 
these  unity  elements  above  the  super  diagonal  is  equal  to  the 
total  number  of  rows  in  row  sets  1  through  (K-l)  and  is  given  by 

K-l 

U  =  E  r (kl)  .  (2.29) 

kl=l 

Substituting  the  value  of  r(kl)  from  (2.28),  and  simplifying, 
equation  (2.29)  becomes 


U 


K-l 
E  (N, 
kl=l 


K-kl+1 


{ 


K 

n  (N.+i)}] 

i=K-kl+2 


(2.30) 
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It  can  also  be  noted  that  the  number  of  columns  in  column  set 


k2,  (k2  =  l,2f...,K)  in  (I-P) ,  is  equal  to  the  number  of  rows 

in  row  set  (K-k2+l)  which  is  given  by  r(K-k2+l)  in  equation 

th 

(2.28).  Therefore,  the  unity  element  in  the  j  column  of  column 

set  k2,  (k2  ®  2,3,...,K),  above  the  super  diagonal  of  (I-P)  is 
th 

located  in  the  j  row  of  row  set  (K-k2+l) .  The  super  diagonal 
of  (I-P)  now  contains  all  non-zero  elements. 

2. 4. 2. 3  Arrangement  of  the  States  of  (I-P)^ 

Tl.e  matrix  (I-P)^  is  obtained  from  (I-P)  by  deleting  the 
first  column  and  the  last  row  of  (I-P).  The  first  column  of 
(I-P)  is 


Vi.  ■  «r1>N2-N3 . V 


and  the  last  row  is 


xJJ*1  =  (0,0,.  ...0) 

These  are  deleted  and  the  resulting  (I-P)^  matrix  has  most  of 
its  elements  in  a  left  triangular  structure  with  non-zero  ele¬ 
ments  along  the  main  diagonal  and  some  unity  elements  above  the 
main  diagonal. 
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Because  Che  lasc  row  and  Che  firsc  column  are  removed  from  Che 
row  see  K  and  Che  column  seC  1  respecCively  of  (I-P)  Co  gee 
(I-P)*,  (2.28)  is  revised  Co  give  Che  number  of  rows  in  row  seC 
kl  of  (I-P)*  as  follows: 


f 


NK-kl+l{ 


K 

n  (n  +i)} 

Jl=K-kl+2 


r' (kl)  =  l 


if  kl  -  1,2, . . . ,K-1  (2 


K 

N  (  n  (N.+l) }-l  if  kl  =  K. 
1  1=2 


The  number  of  columns  in  column  sec  k2  of  (I-P)*  is  equal  Co 

Che  number  of  rows  in  row  seC  (K-k2+l)  of  (I-P)*,  i.e.,  r'(K-k2+l). 

Now,  wich  chis  arrangemenC  of  row  and  column  scaCes  in  (I-P)*" 

and  (I-P),  Che  sCaCes  corresponding  Co  Che  elemenCs  of  Che  sCeady 

scace  deparCure  probabiliCy  veccor  A  in  equacion  (2.19)  are  in 

Che  same  order  as  chac  of  Che  row  sCaCes  of  (I-P)*  and  (I-P)*  Thac 

is,  elemenC  ,  j  =  l,2,...,m,  corresponds  Co  Che  row  sCaCe 

of  (I-P)*  and  (I-P)  and  A  . ,  corresponds  Co  X®**  =  (0,0 . 0) 

m+l  — n 

of  (I-P). 


2.4.3  Inversion  of  a  Lefc  Triangular  Macrix 

The  lefc  Criangular  pare  of  (I-P)*  macrix  can  be  inverCed 
using  Che  expliciC  recursive  expressions  developed  by  SchwarCz 


31) 
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et  al.  [SCHR  73]  for  the  elements  of  the  inverse  of  a  left  tri¬ 
angular  matrix.  Let  C  be  a  left  triangular  matrix  of  size  mxm 

with  elements  c.  . ,  (i,j  -  l,2,...,m),  such  that  c,  .  j  0  for 
1  *  J  J  >  J 

all  j .  If  matrix  G  is  the  inverse  of  C ,  then  the  elements  of  G 
are  given  by,  for  i,j  =  l,2,...,m  [SCHR  73], 


1/Cj  ,j  *  ifl=j 


’i,j 


-  < 


1/c  [-  E  c  g  ]  ,  if  j+1  <  i  <  m  (2. 

J>J  i»11 


n=j+l 


0  ,  Otherwise  . 


It  can  be  seen  that  G  is  also  a  left  triangular  matrix.  An 

examination  of  (2.32)  makes  it  clear  that  g.  's  can  be  calcu- 

^  >  J 

lated  recursively  starting  with  column  m  of  C  and  proceeding 

left  towards  column  1.  Within  each  column  j,  g  's  are  caleu- 

*  J 

lated  starting  with  i  =  j  and  proceeding  down  towards  i  =  m. 

It  is  necessary  at  this  stage  to  obtain  a  relation  with 

respect  to  g  's  and  certain  elements  of  C,  which  is  very  use- 
*  f  J 

ful  in  the  calculation  of  steady  state  departure  average  prob¬ 
abilities,  to  be  considered  later.  Let  the  last  row  of  C  be 
represented  by  c^  and  the  column  (j  =  1,2,..., in)  of  G  be 

represented  by  g , .  Then  the  product  of  the  vectors  c  and  g. 

J  j 

is  equal  to  zero  for  j  *  1,2, . . . ,m-l,  as  C  •  G  =  1.  That  is 
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m 

E  c  g 
n-i  ™.n6n,j 


0  ,  for  j  =  1,2 . m-1 . 


Rearranging  the  terms  and  simplifying. 


m-1 


.  =  -  l  1  c  g  .  Lfor  1  =  1,2,..., m-1  . 

i.j  c  1  m,n6n, j  *  J 

m,m  n-l 


(2 


2.4.4  Modification  of  the  Inverse 


2.4.4. 1  Preliminaries 

In  this  section  the  inverse  of  the  lower  triangular  part 
of  (I-P)^  obtained  through  the  relations  developed  in  section 
2.4.3  is  modified  to  take  into  account  the  unity  elements  which 
lie  above  the  main  diagonal.  This  modification  procedure  is 
based  on  the  Product  Form  method  of  obtaining  an  Inverse  [HADL  61]. 
In  this  method,  in  general,  if  a  matrix  is  formed  by  replacing 
the  sth  column  c  with  a  in  a  matrix  C  of  size  mxm,  whose  inverse 
C  *  is  known,  then  can  be  obtained  by  performing  the  follow¬ 

ing  steps: 

(i)  Compute  the  column  vector,  ■  C 
(ii)  Form  the  column  vector, 


x  = 


where  y  4  0.  (If  y  -  0,  then  C,  has  no  inverse.) 
s  s  x 


.33) 
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(iii) 

(iv) 


This  method  is  chosen  because  of  its  simplicity  and  the 
potential  for  obtaining  recursive  relations.  The  above  steps 
will  be  referred  to  as  Fundamental  Steps  (i),  (ii),  (iii),  and 
(iv)  in  subsequent  discussions. 

After  the  columns  and  rows  of  the  matrix  (I-P)^  are  grouped 
into  sets  as  per  the  arrangement  given  in  section  2.4.2,  the 
unity  elements  above  the  main  diagonal  of  (I-P)^  are  located 
in  the  columns  in  the  column  sets  2  through  K  with  one  unity 
element  in  each  column.  The  modification  of  the  inverse  is  done 
by  considering  each  column  containing  a  unity  element  above  the 
main  diagonal  at  a  time,  starting  with  the  first  column  of 
column  set  2  and  ending  up  with  the  last  column  of  column  set  K. 
As  the  unity  elements  within  each  column  set  are  located  along 
a  distinct  and  different  line,  the  modification  procedure  is 
divided  into  (K-l)  phases,  one  phase  for  each  column  set.  Phase 
p,  (p  *  1,2, . . . ,K-1) ,  modifies  the  previously  obtained  inverse 
by  taking  into  consideration  the  unity  elements  above  the  main 
diagonal  in  the  columns  of  column  set  (p+1) .  Within  each  phase 


Replace  the  s  column  of  the  identity  matrix  I 
of  size  mxm  with  the  vector  x  to  obtain  a  matrix 
E. 

Obtain  the  inverse  of  as 
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there  are  a  certain  number  of  iterations,  each  iteration  modi¬ 
fying  the  inverse,  taking  into  consideration  one  column  in  the 
corresponding  column  set.  The  number  of  iterations  in  phase  p 
is  given  by  r' (K-p)  in  equation  (2.31)  because  the  number  of 
columns  in  the  corresponding  column  set  (p+1)  is  equal  to  the 
number  of  rows  in  the  row  set  (K-p)  in  (I-P)^. 

Before  describing  the  modification  procedure,  it  is  neces¬ 
sary  now  to  briefly  describe  the  symbols  used  in  this  section. 

(a)  S,r(i,j)  specifies  the  row  number  in  (I-P)^  of 

the  j*"*1  row  of  row  set  i,  (i  =  1,2 . K) .  j 

takes  on  values  from  1  to  r'(i).  f,r(i,  j)  is 
given  by,  for  i  =  1,2,...,K, 

i-1 

«r(i,j)  =  E  r'(il)  +  j  (2.34) 

11=1 

where  r'(il)  is  given  in  (2.31). 

(b)  £c(k,n)  specifies  the  column  number  in  (I-P)^ 

th 

of  the  n  column  of  column  set  k,  (k  =  1,2,...,K). 
n  takes  on  values  from  1  to  r'(K-k+l).  S.c(k,n)  is 
given  by,  for  k  =  1,2,...,K, 

k-1 

lc(k,n)  =  E  r '  (K-kl+1)  +  n  .  (2.35) 

kl-1 

fcc(k,F)  represents  the  column  number  in  (I-P)^ 


of  the  last  column  of  column  set  k.  In  other  words, 
Hc(k,F)  =■  Ic  (k,  r 1  (K-k+1) )  with  F  =  r' (K-k+1). 


To  make  the  expressions  simpler,  Hc(k,n)  is 
written  as  Jin  in  cases  where  it  is  known  clearly 
that  the  column  set  under  consideration  is  k. 

Similarly  HF  stands  for  Hc(k,F). 

(c)  The  row  number  in  (I-P)^  of  the  unity  element  above 
the  main  diagonal  in  column  £.c(k,n)  is  denoted  as 
R(Hc(k,n))  which  is  equal  to  Hr(K-k+l,n)  as  per  sec¬ 
tion  2. 4. 2. 2  and  the  symbol  introduced  in  part  (a). 
Whenever  Hc(k,n)  is  written  as  Hn,  R(Hc(k,n))  is 
replaced  by  R(Hn) . 

(d)  The  matrix  containing  only  those  elements  which 
form  the  left  triangular  structure  of  (I-P)^  is 
represented  by  C  and  its  inverse  by  G. 

(e)  The  matrix  which  contains  all  elements  of  C  and  all 
the  unity  elements  above  the  main  diagonal,  from 
left  up  to  and  including  the  one  being  considered 
in  the  ith  iteration  of  phase  p,  is  denoted  as  CP’^ 
and  its  inverse  as  G13'''".  Therefore,  (I-P)^  is  C^K 
as  there  are  a  total  of  (K-l)  phases  and  as  F  here 
stands  for  the  number  of  columns  in  column  set  K 
which  are  taken  into  consideration  in  phase  (K-l). 

At  this  stage  it  is  necessary  to  prove  the  existence  of 

D  •  i 

inverses  for  the  matrices  CK‘  for  1  <_  p  <_  K-l  and  1  <  i  <  r’(K-p). 
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The  following  points  [HADL  61]  are  used  in  this  connection  with 
respect  to  column  vectors. 

(i)  A  set  of  m  column  vectors  a i  *  1,2,..., m,  from 

the  m  dimensional  euclidean  space,  Em,  is  said  to 

be  linearly  independent  if  the  only  set  of  a for 
m 


which  I  a. a.  =0  holds  is  a, 
i-1  i_i  1 


, .  *  a  »  0. 
m 


(ii)  All  the  columns  of  a  nonsingular  matrix  form  a 
linearly  independent  set. 

(iii)  A  set  of  m  linearly  independent  column  vectors  a^» 
i  •  l,2,...,m,  which  spans  Em  forms  a  basis  for  Em. 

(iv)  Any  column  vector  in  Em  can  be  expressed  as  a 

linear  combination  of  the  column  vectors  a^,  i  «=  1,2, 


*nit 


That  is,  b 


m 
i=l 


(v)  Given  a  column  vector  b  f  0  in  Em.  Then,  if  in  the 

expression  of  b  as  a  linear  combination  of  the  basis 

m 


vectors  a.,  i.e.,  b  =  E  6. a. ,  any  vector  a.  for 

i=l  1 

which  (5  +  0  is  removed  from  the  set  a^,^,...,^ 

and  b  is  added  to  the  set,  the  new  collection  of  m 
vectors  is  also  a  basis  for  Em  and  therefore  linearly 
independent. 

Any  left  triangular  matrix  with  all  nonzero  elements  along 
the  main  diagonal  has  an  inverse  [SCHR  73].  With  the  arrangement 
of  the  column  and  row  states  as  per  section  2.4.2,  the  matrix  C 
containing  all  the  elements  of  (I-P)^  which  form  the  left  tri¬ 
angular  structure,  has  non-zero  main  diagonal  elements  and 
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therefore  has  an  inverse.  So  its  m  columns  are  linearly  inde¬ 
pendent  and  form  a  basis  for  Em  as  per  (ii)  and  (iii).  In 
iteration  1  of  phase  1,  the  inverse  is  obtained  for  which 

is  the  same  as  C,  but  with  a  unity  element  above  the  main  di¬ 
agonal  in  column  Hc(2,l),  which  is  column  1  in  column  set  2. 

Let  a  .  represent  the  column  2c (2,1)  of  C.  Because  C  is 

x,C  1 1 ) 

a  left  triangular  matrix,  elements  from  rows  2c(2,l)  to  m  are 

the  only  non-zero  elements  of  a^^  i)  *  Let  k  represent  the 
1:1 

column  He (2,1)  of  C  *  .  Now  b^  has  all  elements  of  i) 

with  an  additional  unity  element  in  row  R(2c(2,l)).  As  per 
(iv),  b  can  be  expressed  as  a  linear  combination  of  the  m 
columns  of  C  and  let  the  corresponding  coefficients  be  denoted 
as  i  =  l,2,...,m.  Because  the  elements  in  rows  2c(2,l)  to 

in  of  b  and  a.^^  i)  are  eclual»  and  because  column  a^^ ^  can¬ 
not  be  represented  as  a  linear  combination  of  the  other  (m-1) 
columns  of  C  as  per  (i) ,  ^)  cannot  be  equal  to  0.  As 

per  (v),  if  a^c^2  is  replaced  by  b  in  the  original  set  of 

m  columns  of  C,  then  the  new  set  Is  linearly  independent. 

1:1 

Therefore  the  matrix  C  ,  which  has  the  new  set  of  m  columns, 
has  an  inverse. 

In  each  iteration  of  each  phase,  a  unity  element  is  added 

to  C  in  different  rows.  Therefore,  following  the  same  arguments 

given  before,  the  existence  of  inverses  GK‘  for  all  p  and  i 

can  be  proved  by  induction.  This  ensures  that  the  values  of 

y  0  in  Fundamental  Step  (ii). 
s 
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Instead  of  going  through  each  iteration  in  each  phase  to 
obtain  the  final  inverse  G^K  Df  (i-p)^,  the  following  ap¬ 

proach  is  used: 

(a)  In  phase  1,  Fundamental  Steps  (i)  through  (iv) 

are  performed  in  iterations  (1)  and  (2)  and  the 

1*1 

expressions  of  the  elements  of  the  inverse,  G 
1:2 

and  G  are  obtained. 

(b)  Based  on  the  results  in  (a),  general  expressions 
for  the  elements  of  the  inverse  G^‘n,  which  is  ob¬ 
tained  after  the  nC^  iteration,  are  developed  and 
proved  by  mathematical  Induction. 

(c)  Using  the  results  of  (b) ,  expressions  are  written 

1 ;  F 

for  G  '  ,  the  inverse  obtained  after  all  iterations 
in  phase  1. 

1  z  F 

(d)  By  comparing  the  elements  of  G  and  G  *  ,  the  ex- 

2  :F 

pressions  for  the  elements  of  G  *  are  developed. 

1:F 

(e)  Based  on  the  expressions  for  the  elements  of  G 

2:F 

and  G  ,  general  expressions  for  the  elements  of 
k:  F 

G  ,  the  inverse  obtained  after  k  phases,  are 
developed.  These  are  proved  by  mathematical  induc¬ 
tion. 

k:F 

(f)  Based  on  the  expressions  for  the  elements  of  G  *  , 

the  expressions  for  the  elements  of  the  final  in- 
verse  G  are  developed. 
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These  steps  are  explained  in  the  following  subsections.  It 

is  helpful  to  recall  now  the  relations  g  =  1/c  .  if  i  ■  j , 

*>J  J  f J 

and  g.  .  »  0  if  i  <  j  from  (2.32).  These  relations  are  used  in 
*»  J 

many  simplifications  in  the  following  subsections. 

2.4.4. 2  Inverses  After  Iterations  (1)  and  (2)  of  Phase  1 

In  phase  1,  the  inverse  G  of  the  matrix  C  is  modified,  tak¬ 
ing  into  consideration  the  unity  elements  above  the  main  diagonal 
in  columns  of  column  set  2.  There  are  r’(K-l)  iterations  in  this 
phase  corresponding  to  the  r'(K-l)  number  of  columns  of  column 
set  2. 

Iteration  (1) : 

In  this  iteration  G  is  modified  because  of  the  unity  element 
above  the  main  diagonal  in  the  first  column  of  column  set  2  which 
is  the  column  £.c(2,l)  of  (I-P)^.  The  unity  element  in  this  column 
is  located  at  position  £r(K-l,l).  So,  with  reference  to  the 
Fundamental  Steps,  £  is  the  column  £c(2,l)  of  (I-P)^  and  s  *  ilc(2,l). 

Fundamental  Step  (i): 


r  Gi 


if  i 


G  is  a  left  triangular  matrix  with  its  elements  g  «  0 

i  >  J 

"  j  from  (2.32).  The  elements  of  £,  which  is  column  £c(2,l) 


of  (I-P)\  are  given  by 

if  i  -  R(U) 

if  i  >  £1  (2.36) 

Otherwise 

where  £1  *  £c(2,l)  and  R(£l)  =  £r(K-l,l).  The  i*1*1  element  of  £ 
is  obtained  from 


m 

yi  =  »  for  1  =  1«2 . m* 


Substituting  the  values  of  qi  from  (2.36)  and  using  the  relation 

g  =  0  if  i  <  j,  the  values  of  y  are  given  by 
if  j  i 


yn-i  i  yn+i  _  ][■* 

y£l  y£l  yU  y£l 


From  (2.37)  and  recalling  that  g 


j>3 


1/c.  .  for  all  j  from  (2.32), 
3  >3 


61 


ym  “  [8u,r(hi)  +  ^ 


The  it^1  element  of  x  is  now  given  by 


[8H1,R(J11)  +  ^ 


if  i  =  ill 


x,  a  ( 


(2.38) 


[81,R(U)  +  iftl8i»JcJ,a] 
[8m,R(<tl)  +  ^ 


,  Otherwise  . 


At  this  stage  it  is  helpful  to  introduce  additional  nota¬ 
tions  to  simplify  the  task  of  representing  different  quantities. 
Equation  set  (2.38)  can  be  rewritten  as 


x^  =  -b(2,Jll,i)  ,  for  i  ■  l,2,...,m  , 


(2.39) 


where 


b(2,U,i) 


[rl(2,U,i)  +  f(2,U,i)l 
v(2,Jtl) 


(2.40) 


where 


rl(2,Jll,i) 


8i,R(U)’  if  1  -  1  <  U 
or 

Hi  <  i  <_  m 
0  ,  if  i  -  HI 


(2.41) 
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f(2,£l,i)  - 


f  -  1  ,  if  i  -  £1 


< 


i 

^  1C1  £1* 
j-£l  J,x,x 


Otherwise 


(2.42) 


and 


v(2,U)  =  ®£l,R(£l)  +  1 


(2.43) 


In  the  notations  just  introduced  the  first  index  within 
the  brackets  signifies  the  fact  that  column  £1,  being  considered 
in  this  iteration,  is  a  member  of  column  set  2. 


Fundamental  Step  (iii) : 

til 

E  is  obtained  by  replacing  the  £ ll  column  of  I  with  x.  The 
elements  of  E  are  now  given  by 


i»  j 


-b(2,£l,i) ,  if  j  =  £1  and 

i  *  1,2,... ,m 

1  ,  if  j  ■  l,2,...,m,  but  +  £1,  (2.44) 

and  i  *  j 

0  ,  Otherwise 


where  b(2,£l,i)  is  given  by  (2.40). 
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Fundamental  Step  (iv) : 


G1:1  =  E  •  G 


where,  as  per  the  notation  introduced,  G^'^  is  the  inverse  of 
C^'^"  which  has  the  same  elements  as  C  and  an  additional  unity 
element  in  column  4c(2,l)  and  row  £r(K-l,l).  The  elements  of 
are  obtained  from 


8.1  j  —  £  ®4  n&n  4  »  if  —  1,2, . . .  ,m 


1:1 


Substituting  the  values  of  e.  „  from  (2.44),  the  values  of  g . ' . 
are  now  given  by,  for  j  =  1,2,..., m. 


1:1 

5i,j 


-ga  jb(2,U,i)  ,  if  i  =  41 


=  < 


8i,j  “  8U  ,  Otherwise 


(2.45) 


Equation  set  (2.45)  can  be  rewritten  as,  for  j  =  l,2,...,m. 


f  -gu>jd(2,U,U,i)  ,  if  i  =  u 


1:1  <• 
g  = 
it  j 


(2.46) 


8i,j  "  8il  j **(2,41,41,1)  ,  Otherwise 
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where 


d(2,n,U,i)  -  b(2,i,l,i)  ,  for  i  -  1,2 . m  (2.47) 

d(2,fcl,Al,i)  stands  for  the  coefficient  of  g(1  ,  in  the 
expression  for  the  element  located  at  the  ith  row  and  the  j C 
column  of  the  inverse  matrix,  obtained  after  including  the 
unity  element  above  the  main  diagonal  of  the  column  £1  of 
(I-P)\  which  is  also  a  member  of  the  column  set  2.  Here  Jll  ■ 

*c(2,l). 

It  may  be  observed  here  that  in  the  calculation  of  the  ele- 
1*  1 

ments  of  G  '  ,  the  structure  of  the  previous  inverse  G  was  used 

to  the  extent  that  g .  .  =  0  if  i  <  j . 

1»J 

Iteration  (2) : 

1:1 

In  this  iteration  G  '  is  modified  considering  the  unity 
element  above  the  main  diagonal  in  column  ic(2,2)  of  (I-P)^" 
which  is  column  2  in  column  set  2.  So  here  is  the  column  12 
of  (I-P)\  s  =  12,  and  R(£,2)  =  £r(K-l,2),  where  12  stands  for 
M2, 2). 

Fundamental  Step  (i) : 


Z 


G 


1:1 


a 


The  elements  of  £  are  given  by 


1  ,  if  i  -  R(£2) 


Ci,£2  *  lf  1  i  *2 


(2.48) 


0  ,  Otherwise 


Yj  *  E  g.'.  q,  ,  i  *  1,2 . m 

J-l  J 


Substituting  the  expressions  for  g  *  and  q,»  and  simplifying 

J  j 

using  12  >  ill  and  noting  that  ^  »  0  if  j  >  £.1,  y^  is  given 


"  R(g, 2)d ^ »  If  i  *  ill 


(2.49) 


gi,R(£,2)  '  d(2,£1,£1,1)gill,R(Jl2)  +  j^2gi,jCj,il2  » 


Otherwise 


Fundamental  Step  (ii): 

The  elements  of  x  are  given  by 


,  if  i  ■  Jl2 


,  Otherwise  • 
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Substituting  the  values  of  from  (2.49)  and  using  the  notation 
in  iteration  (1) , 

■  -b(2,£2,i)  ,  for  i  ”  1,2,. . . ,m  .  (2.50) 

The  index  2  signifies  the  fact  that  column  set  2  is  under 
consideration. 

In  (2.50) 


b(2,£2,i) 


[rl(2,&2,i)  +  f (2,&2,i) ] 
v(2,*2) 


(2.51) 


where 


rl(2,A2,i) 


8i,R(£,2)~  d(2,U,£’1’i)8U,R(i2), 

if  i  *  l,2,...,m  but  *1,  12 
<  -d(2,u,n,i)gUjR(£2)  ,  if  i  -  u 

0  ,  if  i  -  12 

f  -1  ,  if  i  =  £2 


(2.52) 


f(2,*2,i)  -  < 


I  Z  5.2 »  Otherwise 

{  j-112  1,J 


(2.53) 
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and 


v(2,£2)  -  8j12,r(£2)  d(2,U,U,*'2)gU,R(Z2)  +  1  '  (2,54) 

Fundamental  Step  (iii): 

E  is  obtained  by  replacing  column  2,2  of  I  with  x.  The  ele¬ 
ments  of  E  are  given  by 


f  -b(2,2,2,i)  ,  if  j  =  12  and  i  «  1,2,. ...m 


e 


< 


1  , 


if  j  =  1,2, .. .  ,m,  but  i*  £2, 
and  i  =  j 


(2.55) 


0  ,  Otherwise 


where  b(2,£2,i)  is  defined  in  (2.51). 


Fundamental  Step  (iv): 


Here  G1'2  is  the  inverse  of  C1-2,  which  is  similar  to  C1:1, 
but  also  having  a  unity  element  in  column  12  ■  2,c(2,2)  and  row 
R( 8.2)  -  lr(K-l,2). 

1:2 

The  elements  of  G  '  are  given  by 
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for  i,j  -  1,2 . m. 


1:2 


g 


i.j 


m 

Z  e 
2-1 


1,2 


1:1 

S2,j  * 


l- 1 

Substituting  the  values  of  e  and  g  *  , 

2,2  2 ,  j 

(2.46)  respectively,  and  simplifying,  for 


given  by  (2.55)  and 
all  j, 


1:2 

8i,j 


< 


gi,j"  gUJ^(2, 21,21,1)  -  b(2,22,i)  d(2,21,21,22) } 

-  g 12,1 (b(2, 22,i) },  if  1  <  i  <  m 

but  4  21,  22 

_g 2i ,  j  (d ( 2 •  21 ,  21 ,  21)  -  b(2, 22,  21)  d(2, 21,  21,  22)  } 

~822, j {b(2»22,2!)},  if  i  -  21 

-g)lljj{-b(2, 22,22)  d(2,21, 21,22)) 

-  g£2>j{b(2,22,22)},  if  i  -  22 


(2.56) 


where  21  -  2c(2,l)  and  22  =  2c(2,2).  Equation  set  (2.56)  can  be 
rewritten  as,  for  j  -  l,2,...,m, 


,  2c(2,2) 

f  -  I  g„  ,  d(2,2,22,i)  ,  if  i  -  2c(2,l),  2c(2,2) 
2-2c(2,l) 


2c(2, 2) 

g,  ,  -  Z  g p  ,  d(2,2,22,i),  Otherwise 
1,J  2-2c(2,l) 


(2.57) 
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where 


-b(2,£2,£2)  d(2,  SI,  £1,  £2)  ,  if  i  «*  £2 


<1(2, U, £2, i)  =  < 


(2.58) 


d(2,£l,£l,i)  -  b(2,£2,i)  d(2,£l,£l,£2) ,  Otherwise 


and 


d(2,£2,£2,i)  =  b(2, £2,i) ,  if  i  =  1,2 . m.  (2.59) 

Here,  £1  <  £2. 

2. 4. 4. 3  Inverse  After  the  Final  Iteration  of  Phase  1 

By  examining  the  expressions  for  the  elements  of  the  inverse 
and  the  different  quantities  obtained  after  Iterations  1  and  2, 
it  is  possible  to  develop  expressions  for  these  quantities  after 
the  final  iteration  r'(K-l)  in  phase  1.  Mathematical  induction 
will  be  utilized  for  this  purpose  with  the  following  conjecture. 

Conjecture  2.1:  The  elements  of  the  inverse  G^’n,  obtained  after 
the  nth  iteration  (n  <  r'(K-l))  in  phase  1,  are  given  by,  for 
j  =  1,2, .. . ,m, 
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rl(2,An,i)  -< 


Jln-1 

-  lmln  8 


A,R(An) 


d(2,A,An-l,i), 


if  1  <  i  <  A1  or  An  <  i  <  m 


An-1 
-  E 
£=£1 


(2.63) 


gA,R(An) 


d(2, A, Jln-1, i)  ,  if  A1  <  i  <  An 


0  ,  if  i  =  An 


f (2, An,i)  =  < 


-1  ,  if  i  =  An 


2  ,c.  ,  Otherwise 

j-An  1,3  3’  n 


(2.64) 


and 


An-1 

v(2.  An)  =  ginjR()ln)  "  {Jl=^18A>R(An)  d<2>2'>,l«-1.*'n)  J  +  1  *  (2.65) 

In  these  expressions,  An  =  Ac(2,n). 

Now  it  will  be  shown  that  the  expressions  given  in  Conjecture 

s  t 

2.1  are  true  for  the  (n+1)  iteration  if  they  are  true  for  the 
ntl1  iteration. 

Iteration  (n+1) : 

In  iteration  (n+1) ,  the  unity  element  located  at  row 
Ar(K-l,n+l)  of  column  Ac(2,n+1)  of  (I-P)^  is  considered.  So, 
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here  £  is  column  £c(2,n+l)  of  (I-P)J 


Fundamental  Step  (i) : 


„l:n 

y  =  G  £ 


1  ,  if  i  -  RUn+1) 


ci,£n+l  ’  if  1  -  £n+1 


(2.66) 


0  ,  Otherwise  . 


Using  (2.66)  and  (2.60),  the  elements  of  are  obtained  as 


fc«n8£’R(£n+1)d(2,*,*n,i)’  ±f  U  ”  1  1  lU 


(2.67) 


in  i 

8i,R(in+l)"^^18i,R(«n+l)d(2’£'to’i)  +  j.^+18i, jCj ,  ta+1* 

Otherwise  • 


Fundamental  Step  (ii): 


The  elements  of  x  are  obtained  from 
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1 

yln+l 


if  i  =  S.n+1 
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-1  ,  if  i  ■  J.n+1 


f  (2,Jln+l,i) 


(2.70) 


Z  8i  in+1  »  Otherwise 
J-tefl  1,3  3 ’  nl 


v(2,«,n+l)  8iln+ijR(ji,n+i)"^=^i8i>R(i!,n+l)d(2'i!,,in,An+1)}+1  *  (2,71) 


Fundamental  Step  (ill) : 

E  is  formed  such  that 


-b(2,£n+l,i)  ,  if  j  -  £n+l  and  i  -  l,2,...,m 


>i  j  =  <  1  ,  if  j  =  1,2 . m  but  1*  £n+l  and  i  -  j  (2.72) 


0  ,  Otherwise  , 


where  b(2,fcn+l,i)  is  given  in  (2.68). 


Fundamental  Step  (iv) : 


G1;n+1  «  E  .  G1:n 
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Substituting  the  values  of  g}"?  and  e.  from  (2.60)  and 
(2.72),  respectively,  and  rearranging  the  terms,  the  elements 
of  are  given  by,  for  j  *  1,2,. ..,m. 


n+1) 

g  d(2,£,£n+l,i)  ,  if  £c(2,l)  <  i  <.  £c(2,n+l) 

2,1) 


(2.73) 

ic(2,n-vl) 

Z  g .  ,d(2,£,£n+l,i) ,  Otherwise 
t=SU:(2,l)  ,3 


where 


Comparing  equation  sets  (2.68)  through  (2.74)  with  the  correspond¬ 
ing  sets  given  in  Conjecture  2.1,  it  is  clear  that  the  conjecture 

s  t 

holds  for  the  elements  of  the  inverse  after  (n+1)  iteration, 

£h 

if  it  is  true  for  n  iteration.  From  the  results  obtained  after 
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Iterations  1  and  2,  it  is  obvious  that  the  conjecture  is  true 
if  n  ■  1,2.  Hence,  by  induction  it  is  true  for  n  ■  1,2,...,F 
where  F  -  r'(K-l). 

Now  using  Conjecture  2.1,  the  expressions  for  the  elements 
1 J  F 

of  the  inverse  G  ,  obtained  after  all  the  iterations  of  phase 

1,  are  developed.  At  this  stage  it  is  useful  to  introduce  an¬ 
other  quanitty  dl(k,£,i)  which  is  given  by 

dl(k,£,i)  -  d(k,£,£F,i)  (2.75) 

1*F 

The  elements  of  G  can  be  written  as,  for  all  j  «  1,2 . . 

£c(2,F) 

Z  g.  dl(2,£,i)  ,  if  £c(2,l)  <  i<  £c(2,F) 

£«£c(2,l) 

(2.76) 

£c(2,F) 

g  -  Z  g  dl(2,£,i)  ,  Otherwise 

1,3  £*>£c(2,l) 

where  dl(2,£,i)  is  given  by  (2.75)  and  can  be  recursively  found 
using  equation  sets  (2.61)  through  (2.65). 


77 


2. 4. 4. 4  Inverse  After  All  Iterations  of  Phase  2 


1:F 

In  this  phase,  G  ,  the  inverse  obtained  in  phase  1,  is 
modified  to  take  into  account  the  unity  elements  above  the  main 
diagonal  in  the  columns  of  column  set  3.  There  are  r'(K-2)  num¬ 
ber  of  iterations  in  this  phase,  one  each  for  each  column.  In 
each  iteration  all  the  four  Fundamental  Steps  are  carried  out. 

Instead  of  going  through  all  iterations,  it  is  possible  to  de- 

2:F 

rive  the  expressions  for  the  elements  of  G  ,  the  final  inverse 
at  the  end  of  phase  2,  by  noting  the  following  points. 

(i)  The  same  type  of  operations  are  performed  in  each 
iteration  of  phase  2  as  in  the  case  of  phase  1. 

(ii)  In  deriving  the  expressions  of  the  inverse  in 

each  iteration  of  phase  1,  the  structure  of  the 
inverse  matrix  obtained  before  is  utilized  to  the 

extent  that  g .  =  0  if  i  <  j. 

J 

(ill)  The  structure  of  the  column  £  is  utilized  in 

each  iteration.  The  pattern  of  the  elements  in 
columns  of  column  set  3  is  similar  to  that  of 
column  set  2.  Only  the  position  of  the  unity 
element  and  the  first  non-zero  element  after 
unity  is  different  in  each  column,  which  is  re¬ 
flected  in  the  expressions  of  the  inverse  ele¬ 
ments  . 


Now  based  on  these  points  and  equation  set  (2.76),  the  ele- 

2:F  1:F 

ments  of  G  can  be  written  in  terms  of  the  elements  of  G 


as,  for  j  *  1,2, ... ,m, 


Ae(3,F) 

L  g*,  dl(3,£,i) ,  if  £c(3,l)  <  i  <  £c(3,F) 
£=£c(3,l)  *,:1 


(2.77) 

1  *  v  lc(3,F)  . t 

g/.  -  2  g,',  dl(3,£,i)  ,  Otherwise 

1,3  £-£c(3,l)  ,3 


where  dl(3,£,i)  **  d(3,£,!!,F,i) . 

ItF  liF 

After  substituting  the  values  of  g. *  and  g  *  for 
£c(3,l)  <_  £  £  £c(3,F)  from  (2.76)  and  rearranging  the  terms, 
equation  s^t  (2.77)  becomes,  for  all  j, 


2:F 


g 


t.j 


£c(2,F)  lc(3  F) 

g.  .  -  E  g,  ,{dl(2,£,i)  -  l  dl(3,£',i)dl(2, £,£’)} 

1,3  £«£c(2,l)  £,»£c( 3,1) 


£c(3,F) 

l  g  .di(3,£,i) 
£=£c(3,l)  *’»3 


if  1  <  i  <  £c(2,l) 
or  £c(3,F)  <  i  <  m 


£c(2,F)  £c(3,F) 

2  g  {dl(2,£,i)  -  Z  dl(3,£' ,i)dl(2,£,£r)  } 
£=£c(2,1)  X"3  £'-£c(3,l) 

£c(3,F) 

-  Z  g  dl(3,£,i),if  £c(2,l)  <  i  <£c(2,F) 
£-£c(3,1) 

(2.78) 


£c(2,F)  £c(3,F) 

Z  g.  Z  dl(3,£’,i)dl(2,£,£')  } 

£-£c(2,l)  x'’3  £’-£c(3,l) 

£c(3,F) 

Z  g.  dl(3,£,i)  ,if  £c(3,l)  <  i  <  £c(3,F)  . 
£-£c(3,1) 
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In  a  similar  way,  the  expressions  for  d(3,i,in,i),  b(3,in,i) 
and  other  related  quantities  can  be  obtained  on  the  basis  of 
equation  sets  (2.61)  through  (2.65).  In  these  expressions  in 
stands  for  ic(3,n).  The  expressions  are 


f  0  ,  if  i  >  in  and  i  =  1,2 . m 


d(3,i,in,i) 


b(3,in,i),  if  i  =  in  and  i  m  1,2 . m 

(2.79) 

-b (3, in, in)  d(3, i, in-1, in) ,  if  i  <  in 

and  i  *=  in 


[  d(3,i, in-1, i)  -  b(3,in,i)d(3,i,in-l,in)  , 

Otherwise 


where 


b(3,in,i) 


v'(37inT  +  f(3,2n,i)]  , 


(2.80) 


where 


1:F 

8i,R(£n) 


1  •  F 

i=a  8£'R(in)  d(3>Jl'to"1,i)’ 

if  1  <  i  <  ic(3,l)  or  in  <  i  <  m 


rl(3,in,i)  =  \ 


in-1  .  _ 

i=il  d(3'^',£n"1'i), 

if  ic(3,l)  <  i  <  in 


(2.81) 


0 


if  i  ■  in 
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-1 


if  i  *  In 


f(3,£n,i)  =  < 


(2.82) 


j 


i 

E  8 

=tn 


1:F 


cj,£n 


Otherwise 


and 


i • p  in—1  i . p 

''<3't">  -  4»,R<tn)  '  {tfu  4;R«ln)d<3'‘!'£“-1''l”>J+1  •  <2-83> 

Substituting  the  value  of  g^'F  in  (2.81)  through  (2.83),  the  ex¬ 
pressions  become 


£c(2,F) 


^  8o  (dl(2,£,i) 


31*R(in)  £=*c(2,l)  *’R(*n) 


Sln-1 

Z  d(3, £’,£n-l,i)dl(2, £,£')}  - 
£'=£c(3,l) 


£n-l 

2=Ac(3,l)8s,,R(!lTl)  or  in  <  i  <  n 


d(3,£,£n-l,i) ,  if  1  <  i  <  £c(2,l) 


rl(3,£n,i)  =  < 


£c(2,F) 

S.=£c(2,l)8j,,R(Jln) 


{dl(2,£,i)  - 


(2.84) 


£n-l 

Z  d(3,£',£n-l,i)  dl(2, £,£')}  - 

i’-icU.D 

£n-l 

Z  g.  .  .d(3,£,*n-l,i),  if  £c(2,l)  <  i  <  £c(3,l) 

£«2.c(3,l)  *”RUn; 


(Equation  (2.84)  continued  on  next  page.) 
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tc(2,F) 


£n-l 


1  85  p f 0 n\  <  “  z  d(3,£\£n-l,i)dl(2, £,£’)} 

£-Ac(2,l)  ’  1  n;  )l'=ic(3,l) 


2.n-l 

S  g 

S.=  £c(3,l) 


£,R(Jtn) 


d(3,£,£n-l,i),  if  £c(3,l)  <  i  <  £n 


0,  if  i  =  In 


-1  ,  if  1  =  In 


f(3,£n,i)  =  < 


i 

£  g 


[  j=*n 


i,jCj,^n  * 


Otherwise 


(2.85) 


and 


v(3’*n)  g*n,R(Jln)  "  t-tc(2ii)8‘.»(*«){dl<2,1,lB) 


£n-l 

I  d(3,*\iln-l,£n)dl(2,it,Jr)  }  - 
fc'=*c(3,l) 


(2.86) 


Jln-1 

I 

i=£c(3,l) 


g£,R(£n) 


d(3f£,£n-l,£n)  +  1  . 


2. 4. 4. 5  Inverse  After  The  Final  Phase  (K-l) 

Instead  of  going  through  all  the  operations  in  all  the 
phases,  the  expressions  for  the  final  inverse  elements  are  ob- 
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Cained  now  using  mathematical  induction.  First,  a  conjecture  is 
made  as  to  the  expressions  for  the  inverse  elements  after  k 
phases  and  then  it  is  shown  to  hold  after  (k+1)  phases.  In  phase 
k,  the  columns  of  column  set  (k+1)  are  taken  into  consideration. 


Conjecture  2.2: 

After  phase  k,  (k  *»  1,2,...,  K-2) ,  the  elements  of  the  in- 
k*F 

verse  G  *  are  given  by,  for  j  =  l,2,...,m, 


k+1  Hc(k+3-kl,F) 

-  Z  {  E  g.  ,h(k+l,k+3-kl,H,i) } , 

kl=2  Jl“£.c(k+3-kl,l) 

if  Hc(2,l)  <  i  <  He (k+1, F) 


k:F 

?i,j 


< 


(2.87) 


k+1  Hc(k+3-kl,F) 

g,  ,  -  Z  (  Z  g.  .h(k+l,k+3-kl,H,i) } , 

kl=2  H=Hc(k+3-kl,l) 

Otherwise 


where 


k+1  Hc(k',F) 

D(k3, £,i)  E  {  E  D(k' , V ,i)T(k3,k' t, l') }, 

k'=k3+l  H,-Hc(k',l) 

if  k3  <  k+1 

h(k+l,k3,H,i)  =<  (2.88) 


0  ,  Otherwise  • 
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rrwymim.  "ig.-- 


The  quantity  h(k+l,k3,£,i)  stands  for  the  coefficient  of 
g  , ,(£c(k3,l)  <_  l  £  £c(k3,F))  in  the  expression  of  the  element 
at  row  i  of  the  inverse  obtained  after  taking  into  consideration 
the  unity  elements  above  the  main  diagonal  in  the  columns  of 
column  sets  2  through  (k+1) . 

In  (2.88),  for  k'  <  k+1, 


dl(k'  ,1  '  ,i)  ,  if  i  £  «.c(k',F) 
or  i  >  «.c(k+l,F) 


D(k’,A\i)  =  < 


(2.89) 


0  ,  Otherwise 


and 


k’-l  £c(k4,F) 

KkS.k'.i,*’)  -  dl(k3,*,Jt’)  I  {  I  dl(k4,U,Jt’) 

k4=k3+l  U=K.c(k4,l) 

(2.90) 

T(k3,k4,£,££) }  . 


As  seen  earlier, 


dl(kV,i)  =  d(k\*\Jl’F,i) 


The  expressions  of  the  other  related  quantities  are  given  in 
the  following  equations.  In  these  equations,  kl  can  take  any 


JNCLASSIFIEO 


SYRACUSE  UNXV  NY  SCHOOL  OF  COMPUTER  AND  INFORMATION  —ETC  F/G  9/2 
MULTIPLE  FINITE  SOURCE  OUEUEING  MODEL  WITH  FIXED  PRIORITY  SCHED— ETC<U) 
JAN  81  M  J  CHANDRA  F30602-77-C-0235 

RADC-TR-80-379-V0L-2  NL 


of  the  values  l,2,...,k+l  and  £n  -  £c(kl,n): 


(i) 


d(kl,£,£n,i) 


where 

(il)  b(kl,£n. 


where 

(ill) 


rl(kl,£n,i)  - 


0,  if  l  >  £n  and  i  »  l,2,...,m 


b(kl,£n,i)  ,  if  £  =  £n  and  i  ■  1,2,. . .  ,m 

<  (2.91) 

-b(kl,£n,£n)d(kl,£,£n-l,£n)  if  l  <  £n 

and  i  ■  Jin 

d  (kl ,  l ,  £n- 1 ,  i  ) -b  (kl ,  in ,  i  )  d  (kl ,  l ,  £n-l ,  in)  , 

Otherwise 


i)  -  7fri7aS)  [rl(kl,£n,i)  +  f (kl,£n,i) ]  ,  (2.92) 


8 


i,R(£n) 


kl  £c(kl+2-il,F) 

E  {  E  8.  »/, n*h'(kl,kl+2-il,£,i»- 

il-3  £=*£c(kl+2-il,l)  x"  ; 


£n-l 


E  8«  tj  /  tn\ ^ (kl , £ , £n-l, i)  , 
£=£c(kl ,  1)  £*RUn> 


if  1  <  i  <  £c(2,l)  or  £n  <  i  <  m 


kl  £c(kl+2-il,F) 

E  (  E  g.  _...h'(kl,kl+2-il,£,i)}- 

il=3  £=£c(kl+2-il,l)  *,KUIW 

£n-l 

E  8  n  Rf  2n-l,i)  , 

£*£c(kl,l) 

if  £c(2,l)  <  i  <  £n 


0  ,  if  i  =  £n 
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where 


(iv)  h'  (kl,k6,2,l)  -  D(k6,i,i)  - 


kl-1  f,c(k',F) 

I  {  Z  D(k'  ,1'  ,i)  T(k6,k’  ,*.,8.’) }-  (2.94) 

k’=k6+l  H'<tc(kM) 


in-1 

1  d(kl,£', An-1,  i)  T(k6,kl,  i  ,£ ' ) , 

«.’=£c(kl,l) 


where  T(k3,k' >£,£')  is  given  by  equation  (2.90)  and  D(k6,£,i) 
is  given  by 
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(viiD  v(kl.An)  -  gln§RUn)  - 


kl  £c(kl+2-il,F) 
E  {  E 

il-3  (kl+2-il ,  1) 


8£,R(in)h'(kl*kl+2-il*4'4n)}- 


(2 


in-1 

1  8i  Rrin\d(kl*£»£n-1,£n)  +  1* 
i-M-(kl.l) 


Now  the  expressions  of  the  elements  of  the  inverse  after 
phase  (k+1)  are  derived  assuming  that  the  expressions  given  in 
Conjecture  2.2  are  true  for  phase  k. 

Phase  k+1: 

In  phase  (k+1) ,  the  columns  of  column  set  (k+2)  are  con¬ 
sidered.  Instead  of  going  through  all  iterations  of  phase  (k+1) , 
the  method  used  to  obtain  the  expressions  in  phase  2  given  the 
expressions  in  phase  1,  can  be  utilized.  The  three  arguments 
given  in  phase  2  are  general  in  nature  which  can  be  applied  to 
any  phase. 

k+l:F 

Now,  based  on  equation  set  (2.76),  the  elements  of  G  , 

the  inverse  after  phase  (k+1),  can  be  written  in  terms  of  the 
k'F 

elements  of  G  as,  for  all  j, 
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£c(k+2,F) 

E  g^:F  dl(k+2,£,i),  if  £c(k+2,l)  <  i  <  £c(k+2,F) 
£=£c(k+2,l)  11,3 


k+1 :  F 

i.j 


g 


(2.98) 


W-F  £c(k+2,F)  , 

g/.  -  E  g  '  dl(k+2,£,i) ,  Otherwise 

£=£c(k+2,l)  ‘3 


where  dl(k+2,£,i)  =  d(k+2,£,£F,i) . 

ki  F 

Now  the  value  of  g . *  can  be  substituted  in  the  preceding 

J 

equations  and  the  resulting  relations  can  be  simplified.  Because 
of  the  lengthy  procedure  involved  in  the  simplification  process, 
the  details  of  simplification  for  the  expressions  in  (2.98)  when 
1  ±  i  <  £c(2,l)  and  £c(2,l)  _<  i  £  £c(k+l,F)  are  given  in  Appendix 
C.  Details  of  simplifications  for  other  ranges  are  not  given  as 
the  procedure  is  similar  to  the  one  given.  Upon  simplification 
and  rearranging  the  terms,  equation  set  (2.98)  becomes 


k+2  £c(k+4-kl,F) 

-  E  (  E  g.  . h (k+2 ,k+4-kl , £ , i) } , 

kl=2  £=£c(k+4-kl,l)  ,3 

if  £c(2,l)  <  i  <  £c(k+2,F) 


k+l:F 

8i,j 


< 


(2.99) 


k+2  £c(k+4-kl,F) 

g,  ,  -  2  {  E  g  h (k+2 ,k+4-kl , £ , i) } , 

kl“2  £=£c(k+4-kl,l)  Jt’3 


Otherwise 
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h(k+2,k3,£,i) 


where 


D(k3, £ , i)  - 


and 

T(k3,k' , £ 


f  D(k3,£,i) 


k+2  £c(k\F) 

Z  {  Z  D(k\£\i)T(k3,k 

k’-k3+l  £'-£c(k',l) 

if  k3  <  k+2 


< 


(2 


0  ,  Otherwise 


/ 

dl(k3,£,i)  ,  if  i  <  £c(k3,F) 
or  i  >  £c(k+2,F) 

< 


(2 


0  ,  Otherwise 


£')  -  dl(k3,£,£')  - 

(2 

k'-l  £c(k4, F) 

Z  {  Z  dl(k4, ££,£') T(k3,k4,£,££) } , 

k4“k3+l  ££-£c(k4,l) 


.100) 


.101) 


.102) 


when  k’ 


k+2. 


The  expressions  for  D(k3,Jl,i)  for  k3  £  k+1  and  T(k3,k'  ,1,1') 


(iii) 


k:F 

8i,R(ta) 


*»-l  k-F 


rl(k+2,^n,i) 


if  1  <  i  <  Jlc(k+2,1) 
or  Jin  <  i  <  m 


- 


Jln-1 


k:F 


if  £c(k+2,l)  £  i  <  Jin 


0  ,  if  i  -  Jin 


(2 


(iv) 


f  (k+2,Jln,i)  -  < 


-1  ,  if  i  ■  Jin 


^  8i|j  Cj,Jln  ’  °therwiSe 
j=£n 


(2. 


and 


Jln-1 


(v)  v(k+2, in)  -  g^R((n)  -  {  I 

1  (2 


k;F 

The  expressions  for  g, '  given  in  (2.87)  can  be  substituted 

i  >  J 

in  the  expressions  (2.105)  through  (2.107)  and  simplified.  The 

same  procedure  is  followed  in  the  derivation  of  the  elements 
k+1  •  F 

g  '  and  so  the  details  are  omitted  here.  Upon  reduction,  the 
i»  J 

expressions  become 


.105) 


106) 


107) 


(ill) 


f 


k+2  £c(k+4-il,F) 

r  {  1  g 

11*3  £■*  ic(k+4-il,l) 


h'(k+2,k+4-il,£,i)} 


g 


i,R(£n) 


£,R(£n) 


£n-l 


£  8a  tj  /  (k+2 ,  £,  £n— 1 ,i)  , 

£=  £c(k+2,l)  ,R(  n; 


if  1  £  i  <  £c(2,l)  or  in  <  i  £  m 


k+2  £c(k+4-il,F) 

rl(k+2,£n,i)  =  <  -  I{  Z  g.  ...  .h’ (k+2,k+4-il,i,i)  (2.108) 


il-3  £=£c(k+4-il,l) 
in-1 

i=£c(k+2,l)  8ji*R^n> 


if  £c(2 ,1)  £  i  <  in 


d (k+2, i, in-1, i) , 


0  ,  if  i  =  in  , 


where 


k+1  ic(k’.F) 

h'(k+2,k6,i,i)  *  D(k6,i,i)  Z  {  Z  D(k* ,i ' ,i)T(k6,k’ ,i,i') 

k ' =k6+l  £,»ic(k',l) 

£c(k+2,F) 

Z  d(k+2, £', in-1, i)T(k6, k+2, £,£')  ,  (2.109) 

£'=£c(k+2,l) 


where  d(k+2, £', in-1, i)  is  given  by  (2.103).  Expressions  for 
D(k6,£,i)  are  given  by  (2.101)  when  k6  =  k+2  and  by  (2.89)  when 
k6  £  k+1.  Expressions  for  T(k6,k' ,’’.,£’)  are  given  by  (2.102) 
when  k’  =  k+2  and  by  (2.90)  when  k'  £  k+1. 
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(iv) 


-1 


if  i  -  in 


f  (k+2,  Jin,  i)  ■  < 


(2.110) 


i 

I 

j*=£n 


8i,jCj.in 


Otherwise 


and 


(v)  v (k+2, in)  =  g£nfR(£n)  - 

k+2  ic(k+4-il,F) 

S  {  S  g0  h’ (k+2,k+4-il,i,in) }  -  (2.111) 

il=3  i=Jlc(k+4-il,l)  ’  1  ' 


in-1 

^  8o  o ( o„^ ^ (k+2 , i , in— 1 , in)  +1. 
i=ic(k+2,l)  *RU  ' 


Comparison  of  the  expressions  obtained  after  phase  (k+1) 
with  the  corresponding  ones  in  Conjecture  2.2  reveals  that  the 
expressions  given  in  the  Conjecture  hold  for  k  “  k+1.  An  exami¬ 
nation  of  the  expressions  obtained  after  phases  1  and  2  makes  it 
clear  that  the  expressions  given  in  Conjecture  2.2  hold  for 
k  =  1,2.  By  mathematical  induction,  then,  Conjecture  2.2  holds 
for  k  -  1,2, . . . , (K-l) ,  where  (K-l)  is  the  total  number  of  phases. 

It  is  now  possible  to  write  the  expressions  of  the  elements 
of  final  inverse  Gv  ",  which  is  obtained  after  modifying  the 
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original  inverse  G  of  the  lower  triangular  part  of  (I-P)\  taking 

into  consideration  all  the  unity  elements  above  the  main  diagonal 
1  k-1 • F 

of  (I-P)  .  The  elements  of  G  ’  are  given  by,  for  j  “  l,2,...,m. 


K  £c(K+2-il,F) 

-  St  S  g„  .h (K ,K+2-il ,  £ ,  i)  }  , 

il“2  £=£c(K+2-il,l) 

if  £c( 2,1)  <  i  <  m 


(K-1) :F 

8i,j 


=  < 


(2.112) 


K  £c(K+2-il,F) 

8,  ,  "  S  {  S  g.  ,h(K,K+2-il,£,i)}, 

,J  il-2  £**£c(K+2-il,l) 


Otherwise 


where,  for  kl  =  2,3,...,K, 


K  £c(k\F) 

(i)  )»lk,kl,£,i)  =  D(kl, £,i)  Z  {  Z  D(k' ,£' ,i)T(kl,k’ ,£,£’) } 

k ' =kl+l  £'=£c(k',l) 

(2.113) 


(ii) 


D(kl,£,i)  =<^ 


dl(kl, £ , i)  ,  if  i  <  £c(kl,F) 


0  ,  Otherwise 


(2.114) 


(iii) 


dl(kl, £,i)  =  d(kl,£,£F,i) 


(2.115) 


The  values  of  all  other  related  quantities  are  given  as 
per  expressions  (2.90)  through  (2.97)  in  Conjecture  2.2  for 
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kl  -  2,3 . K. 

2.5  CALCULATION  OF  STEADY  STATE  PROBABILITIES 

In  this  section  the  elements  of  the  steady  state  departure 
probability  vector  A  are  calculated  using  the  equations  (2.20)  and 
(2.18)  and  the  expressions  of  the  elements  of  the  inverse  of  (I-P)^ 
obtained  in  the  previous  section.  Rewriting  equation  (2.20)  as 

(A1,A2,...,Ajn)  -  Arofl(Z1,E2 . zm)c(  )  ’ 

each  element  of  the  row  vector  on  the  left  hand  side  can  be 
expressed  as 


m 

A  z 
nrfl  i-1  1 


2,  8 


(K-l)  :F 


for  j  =  1,2, . . . ,m  .  (2 


( K— 1 }  •  F 

After  substituting  the  values  of  g  ‘  from  (2.112)  and  re- 

i>  J 

arranging  the  terms,  equation  (2.116)  can  be  written  as 


Aj  ■  Vi(sj  -  SJ> 


for  j  -  1,2, . . . ,ra, 


(2 


where 


x  M2,l)-1 

4  =  z  Z  g  ,  for  j  ■*  1,2 . ra  (2 

J  i=1  1  J-.J 


.116) 


.117) 


.118) 
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K  £c(K+2-il,F)  m 

=  l  [  E  g  1  I  Z  h(K,K+2-il,M)}] 

J  il=2  «.-Hc(K+2-il,l)  ’ 3  i=l 


for  j  «=  1,2, . . . ,m  . 


1  2 

The  expressions  of  and  are  now  modified  to  reduce 
them  as  functions  of  the  quantities  corresponding  only  to  the 
rows  lc(2,l)  to  m  of  (I-P)^".  The  number  of  rows  from  ii,c(2,l) 
to  m  is  less  compared  to  the  remaining  rows  in  most  cases,  and 
this  modification  reduces  the  computer  storage  and  the  amount 
of  computation  time  in  such  cases. 


2.5.1  Calculation  of  S. 


It  is  necessary  at  this  stage  to  consider  a  lower  triangu¬ 
lar  matrix  L  of  size  (m+l)x(ra+l)  whose  elements  are  given  by 


c,  .,  for  i  =  1,2,. ..,m  and  j  ■  1,2,. . . ,m 
i»  J 


(2.119) 


-Z  ,  for  i  =  m+1  and  j  =  1,2 . m 


(2.120) 


1  ,  for  i  =  nrt-1  and  j  =  m+1 


0  ,  Otherwise 
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where  the  c.  ,'s  are  the  elements  of  the  lower  triangular  part 

1  9  J 

1  of 

of  (I-P)  and  the  -Z^ 's  are  the  last  m  elements  of  the  (nrt-1) 
row  of  (I-P). 

Because  of  the  nature  of  the  relationship  between  L  and 

(I-P)1  matrices,  the  values  of  g.  's  related  to  (I-P)1  and  L 

1  >  J 

are  the  same  for  the  corresponding  values  of  i,j  =  l,2,...,m. 
Therefore,  using  (2.33),  the  following  relation  can  be  estab¬ 
lished  with  respect  to  the  elements  of  L: 


nrfl.j 


Z/mfl,i8i,j  f0r  J  *  i.2, . . .  ,m 


(2.121) 


Substituting  the  value  of  ^  from  (2.120)  into  (2.121), 


=  “i8i,j  ’  f°r  j  3 


£c(2,l)-l  m 

1  ^i®i  i  +  Z  Zi8i  i 
i=l  1  1,3  i=ic(2,l)  ’3 


=  St  +  I  Z  g 
3  i*=£c(2,l)  1  1,3 


Therefore 


1 

S,  =  8nr+0  1  "  1  Zl8i  1  * 

j  nr+-l.j  i=e,c(2,l)  1  1,3 


(2.122) 


for  j  »  1,2, ... ,m 
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where 


m 

H(K+2-il, S,)  =*  l  Z  h(K,K+2-il,Z,i)  .  (2.124) 

i=l 


If  (K+2-il)  is  written  as  kl  for  the  sake  of  simplicity,  then 
(2.124)  becomes 


ra 

H(kl,£)  **  l  Z  h(K,kl,i,i)  . 
i-1 


Substituting  the  value  of  h(K,kl,Jl,i)  from  (2.113),  H(kl,H)  is 
given  by 
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m 

H(kl,I)  -  l  Z  (D(kl,*,i)  - 
i-1  1 


K  Jtc(k\F) 

Z  {  r  D(k\*\i)  T(kl,k\ &,£*)}] 

k’-kl+l  £'-Jlc(k\l) 


(2.125) 


-  Z  Z  D(kl,i,i)  - 
i-1 

K  Hc(k*,F)  m 

Z  {  Z  T(kl,k\ «.,«•')[  2  Z.D(k\i\i)]). 

k'=kl+l  Z'-ic(kM)  i-1 


In  (2.125) 


m  Hc(kl,F)  m 

Z  Z  D(kl,£,i)  =  Z  Z  D(kl,l,l)  +  Z  ZD(kl,l,i) 

i-1  1  i=l  i-Hc(kl,F)+l  1 


£c(kl,F) 

Z  Z  dl(kl,fc,i). 
i-1 


because  of  (2.114).  Similarly, 


m  £c(k',F) 

I  Z  D(k',*',i)  =  Z  Z  dl(k',£',i)  . 
i-1  i-1 


Now  it  is  necessary  to  introduce  a  new  quantity  Y(kl,£,£n) 
which  is  defined  as 
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Y (kl, A, An) 


(2.126) 


An 

=  Z  Z  d(kl,A,An,i)  . 
1=1  1 


Then,  because  dl(k’,A’,i)  -  d(k' , A ' ,A 'F,l)  as  per  (2.115),  equation 
(2.125)  can  be  rewritten  as 


K  Ac(k',F) 

H(kl, A)  =  Y(kl,A, AF)  Z  {  Z  T(kl ,k* , A, A ' ) 

k '  *kl+l  A'-Ac(kM) 

(2.127) 

Y(k,,A,,A'F)}  , 


where  T(kl,k' ,A, A')  is  given  by  (2.90),  AF  =  Ac(kl,F)  and  A'F  - 
Ac(k' ,F) . 

In  (2.127)  the  only  unknown  quantities  are  the  newly  intro¬ 
duced  Y(.,.,.)'s,  which  were  defined  in  (2.126).  Substituting 
the  values  of  d(kl,A,An,i)  from  (2.91)  into  (2.126)  and  simplify¬ 
ing,  the  values  of  Y(kl,A,An)  can  be  expressed  as 

0  ,  if  A  >  An 

An 

Y(kl,A,An)  *  ^  Z  Z  b(kl,An,i)  ,  if  A  -  An  (2.128) 

i=l  1 

Y (kl, A, An-1)  -  d (kl , A , An-1 , An) Y (kl , An , An) , 
otherwise. 
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The  expression  set  (2.128)  forms  a  set  of  recursive  relations. 

For  any  value  of  kl,  Jin  can  be  varied  from  £1  to  £F.  For  each 
value  of  £n,  £  can  be  varied  starting  from  £n  and  ending  up  with 
£1.  When  £  ■  £n,  the  value  of  Y(kl,£n,£n)  is  obtained  and  then 
the  values  of  Y(kl,£,£n)  for  £1  <^  £  <  £n  are  determined  recur¬ 
sively  using  (2.128).  So  the  next  step  now  is  to  calculate 
Y(kl,£n,£n). 

Using  the  value  of  b(kl,£n,i)  from  (2.92)  in  (2.128), 

Y(kl,£n,£n)  is  given  by 

Y (kl, £n, £n)  =  ^(klllrO  ^(kl.in)  +  t2(kl,£n)]  ,  (2.129) 


where 


t^(kl,£n) 


£n 

E  Z  rl(kl,£n,i) 
i=l  1 


(2.130) 


and 


a 

t  (kl,£n)  -  E  Z  f(kl,£n,i)  .  (2.131) 

i=l  1 

Now  the  expressions  (2.130)  and  (2.131)  are  to  be  modified  so 
1  2 

that  t  (kl,£n)  and  t  (kl,£n)  can  be  expressed  as  functions  of 
the  quantities  related  to  the  rows  from  £c(2,l)  to  m  only. 

First  t^(kl,£n)  is  considered.  Substituting  the  value  of 
r(kl,£n,i)  from  (2.93)  and  rearranging  the  summations,  (2.130) 
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can  be  rewritten  as 


102 


(ii)  As  per  (2.132) 


ln-1 

1  Z  d(kl,  l,  ln-l,i)  =  Y(kl,  S.,  to-1)  .  (2 

i=l  1 


(iii)  Using  the  values  of  h’ (kl,kl+2~il,£,i)  from  (2.94) 
and  rearranging  the  summations, 


£.n-l  2.n-l 

Z  Z  h'(.a,kl+2-il,K.,i)  =  Z  Z  D(kl+2-il,fc,i) 
i=l  1  i=l  1 


£n-l  £n-l 

Z  T(kl+2-il,kl,i,*.*){  Z  Z . d (kl ,  £ ' , 5-n-l , i)  } 
£»=£c(kl,l)  i=l  1 


kl-1  £c(k\F)  to-1 

Z  [  Z  T(kl+2-il,k' ,  £,£.'){  Z  Z.D(k' 

k'=kl+2-il+l  H'=ic(k',l)  i=l  1 


Because  of  (2.95) 


£n-l  £c(kl+2-il,F) 

Z  Z .  D(kl+2-il,  £,  i)  =  Z  Z.d(kl+2-il,  SL,  £F,i)  , 

i=l  1  i=l  1 


where  S.F  =  £c(kl+2-il,F) . 
As  per  (2.126), 


£c (kl+2-il,F) 

Z  Z  d(kl+2-il,  £,  £F,i)  =  Y(kl+2-il,£,  JLF) 

i=l 


and  so 


134) 


.135) 
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Jin—  1 

E  Z  D(kl+2-il,£,i)  =  Y(kl+2-il,£,£F) . 
i=l  1 


Similarly 


in-1 

E  Z.D(k,,i',i)  =  Y(kV,£'F)  , 
i=l  1 

where  i'F  =  £c(k',F). 

Now,  based  on  the  modifications  in  (i) ,  (ii) ,  and  (iii), 
(2.132)  can  be  rewritten  as 


m 

t  (kl, in)  =  gnrfl)R(ln)  -  is£c(2jl)Zi8i,R(£n) 


£c(kl,n-l) 

^  R( SLn)^  ^ 

J.=£c(kl,l)  *',RUn) 


kl  £c(kl+2-il,F) 

-  E  [  E  g  {Y(kl+2-il,£,£F) 

il=3  £=£c(kl+2-il,l)  n; 


£n-l 

E  T(kl+2-il,kl , l  ,£  ’  )Y  (kl,  5.  ’  >  £n-l) 
l ,=£c(kl,l) 


kl-1  £c(k\F) 

E  (  l  T (kl+2-il ,k ' , £ , £ ' )Y (k ' , £ 

k’=kl+2-il+l  £'=£c( k',1) 


2 

Second,  the  expression  for  t  (kl,£n)  in  (2.131)  is  modified. 
Substituting  the  values  of  f(kl,£n,i)  from  (2.96)  and  simplifying 


(2.136) 


VF))}]  . 


10A 


where  £F  =  £c(kl+2-il,F) ,  £'F  =  £c(k',F),  and  the  values  of 
v(kl,£n)  and  T(. can  be  obtained  from  (2.97)  and  (2.90), 
respectively. 

Now  using  (2.138)  and  the  recursive  relation  set  (2.128), 
the  values  of  Y(.,.,.)  can  be  obtained.  These,  along  with 
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T(.,.,.,.)'s  from  (2.90),  give  the  values  of  H(.,.)  in  (2.127). 

2 

Then  (2.123)  can  be  used  to  obtain  S^. 


2.5.3  Final  Expressions  for  the  Steady  State  Departure  Probabilities 


1  2 

The  values  of  S^  and  can  be  calculated  using  the  expres¬ 
sions  obtained  in  subsections  2.5.1  and  2.5.2  and  then  these 
can  be  substituted  in  (2.117)  to  obtain  the  value  of  the  steady 
state  departure  probability,  ,  (j  *  l,2,...,m),  in  terms  of 

A  . , .  Then  the  normalizing  condition  (2.18)  can  be  used  to  ob- 
nrt-1 

tain  the  value  of  A  ,  j  ■  1,2, . . .  ,ntfl. 

1  2 

Alternately,  instead  of  computing  S^  and  separately  and 

1  2 

using  these  values  in  (2.117),  the  expressions  of  S^  and  S^  in 
(2.122)  and  (2.123),  respectively,  can  be  substituted  in  (2.117). 

Then  the  combined  expression  can  be  used  for  the  computation  of 
Aj .  In  this  way,  A^  can  be  written  as 

K  £c(K+2-il,F) 

A.  "  A  ..[g..  ,  -  l  {  E  g,  ,<Z  +  H(k+2-il,J0)}]  .  (2.139) 

j  m+1  art-1, J  ilm2  j,=£C(K+2-il,l)  * 


Using  the  normalizing  condition, 


E  A  -  1 
j-1  J 


given  in  (2,18),  the  value  of  A 


can  be  obtained  from 


(2.140) 


^m+1 


[1  + 


m  -1 

2  W  ]  \ 

J-l  3 


where 


W  m  2 
j  8m+l,j 


K  Jlc(K+2-il,F) 

Z  {  Z  g  (Z  +  H(K+2-il,2))  }, 

il-2  i=£c(K+2-il,l)  ’3 


(2.141) 


for  j  -  1,2, . . . ,m  . 


In  (2.141),  H(K+2-il,«,)  is  calculated  from  (2.127),  (2.128),  and 
(2.138).  The  values  of  are  now  given  by,  for  j  «  l,2,...,m, 

Aj  -  .  (2.142) 


2.6  SUMMARY  OF  THE  ALGORITHMS 

The  solution  development  of  the  algorithms  was  explained  in 
sections  2.2  through  2.5.  The  expressions  of  the  required  quanti¬ 
ties  were  also  derived.  In  this  section  the  solution  procedure 
of  the  algorithms  is  summarized  in  a  sequence  necessary  to  com¬ 
pute  the  mean  values  of  the  system  performance  measures. 

The  sequence  of  the  algorithms  for  implementation  is  given 
in  Figure  (2.5).  Now  each  algorithm  in  this  sequence  is  described 
in  steps,  organized  in  a  manner  suitable  for  adaption  as  a  computer 
code  in  the  following  subsections. 
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SEQUENCE  OF  THE  ALGORITHMS 


FIGURE  2.5 


t 


2.6.1  Algorithm  for  Generating  the  States  of  (I-P)^- 

The  details  of  the  procedure  of  this  algorithm  which  were 
discussed  in  section  2.4.2  are  now  combined  and  explained  in 
the  following  steps. 

A.  Row  States:  (Numbered  from  the  top) 

Step  1:  (Initialization).  Set  i  •  1,  kl  *  1,  il  -  1, 

£1  ■  K,  and  £2  =  K-l.  Set  the  ordered  row 

sets  1  through  K  empty.  Set  the  first  row 

state  equal  to  X*  =  (N^.N^ . NK-1,NK-1)‘ 

Set  j£  -  N£  for  £  -  1,2 . K-l  and  jR  -  J^-l. 

Add  the  row  state  X^"  as  the.  first  member  of 
— n 

the  ordered  row  set  1.  Set  il  **  il  +  1  and 
i  *  i  +  1. 

Step  2:  If  j£l  ^  0,  set  j£1  -  j £1~1  and  go  to  Step  4. 
Otherwise  go  to  Step  3. 

Step  3:  Set  *  N£^  and  £1  *  £1-1.  Go  to  Step  2. 
Step  4:  Set  the  iC^  row  state, 


( j  1 » j  2  * 


”V’ 


If  either  £2  =  0  or  the  £^  element  of 
=  Ni2’  80  t0  Step  Otherwise  go  to 
Step  6. 


109 


Step  5:  Add  row  state  as  the  il*'*1  member  of  the 

“H  K 

ordered  row  set  kl.  If  i  =  1  H  (N,+1)}-1, 

k=l 

go  to  Step  7.  Otherwise  set  il  =  il+1, 

i  =  i+1,  and  £.1  =  K.  Go  to  Step  2. 

Step  6:  Set  kl  =  kl+1  and  il  =  1.  Add  row  state  X1 
r  ~n 

as  the  ilC^  member  of  the  ordered  row  set  kl. 

Set  2.2  =  2.2-1,  il  =  il+1,  i  =  i+1,  and  21  =  K. 

Go  to  Step  2. 

Step  7:  Stop. 

B.  Column  States:  (Numbered  from  left) 

Step  1:  (Initialization).  Set  v  =  2,  kl  =  1,  and  vl  =  1. 

Set  ordered  column  sets  1  through  K  empty. 

Step  2:  Set  the  vt!l  member  of  row  set  (K-kl+1)  equal  to 
t  h 

the  vl  member  of  the  ordered  column  set  kl. 

If  v  is  equal  to  the  number  of  rows  in 
row  set  (K-kl+1),  then  go  to  Step  3.  Other¬ 
wise,  go  to  Step  A. 

Step  3:  If  kl  =  K,  go  to  Step  5.  Otherwise  set  v  =  1, 
kl  =  kl+1,  and  vl  =  1.  Go  to  Step  2. 

Step  4:  Set  v  =  v+1,  and  vl  =  vl+1.  Go  to  Step  2. 

Step  5:  Stop. 

After  these  steps,  there  are  (tiH-1)  number  of  row  states  and 

K 

m  number  of  column  states,  where  m  is  given  by  {  II  (N.+l)-2). 

2=1  * 

The  m  column  states  and  the  first  m  row  states  form  the  states 
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of  (I-P)*  matrix.  The  (m+l)St  row  state,  i.e.,  X®*1  -  (0,0,...,  0) , 

“ n 

is  required  for  calculating  the  elements  of  the  row  vector  Z  and 
of  the  matrix  L  defined  in  (2.120). 


2.6.2  Algorithm  for  the  Calculation  of  Steady  State  Probabilities 


As  the  expressions  for  the  elements  of  the  inverse  of  (I-P)* 
were  used  to  obtain  explicit  expressions  for  the  steady  state  de¬ 
parture  probabilities,  the  algorithms  for  inversion  of  (I-P)^  and 
for  the  calculation  of  the  steady  state  probabilities  can  be  com¬ 
bined  together  for  the  purpose  of  implementation.  The  recursive 
expressions  for  calculating  the  inverse  elements  do  not  require 
all  the  elements  of  the  transition  probability  matrix  P  at  one 

time.  The  g  '  s  are  calculated  as  and  when  the  elements  of  P 

j 

are  generated.  Therefore,  for  the  sake  of  implementation,  the 
algorithm  for  generating  the  elements  of  P  can  also  be  combined 
with  the  algorithms  for  inversion  of  (I-P)^  and  for  the  calcula¬ 
tion  of  steady  state  departure  probabilities. 

Step  1:  Consider  the  matrix  L,  whose  elements  are  given 
by  (2.120). 


Decrease  j  by  1  from  m+1  until  it  is  equal 

to  1.  For  each  value  of  j  increase  i  by  1  from 

j  to  m+1.  For  each  value  of  the  pair  (i,j) 

calculate  the  values  of  the  elements  p.  's,  cor- 

J 

responding  to  the  row  and  column  states,  using 

the  formulae  given  in  section  2.3,  and  Jt  's 

1 » J 
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Step  2 
Step  3 


Step  4: 


using  (2.120).  Using  the  recursive  relations  (2.32),  cal¬ 
culate  g.  's  one  by  one  recursively,  with  £.  .  m  c.  .. 

1  #  J  1  >  J 

Store  the  values  of 

(i)  g .  '  s  for  He (2,1)  <_  i  <  ro+1  and  1  <  J  <  nr+1 

l,  j 

(ii)  £  's  for  i  =  1,2 . m  and  j  =  i 

i»  3 

(iii)  £  's  as  c  's  for  £c(2,l)  i,j  £  m,  and 

i» j  i » 1 

(iv)  £.  's  as  -Z  's  for  i  =  m+1  and  1  £  j  <_m. 

i »  j  3 

Set  kl  =  2. 

Increase  £n  by  1  from  £c(kl,l)  to  £c(kl,F).  For 
each  value  of  £n,  (i)  increase  i  by  1  from  £n  to 
m,  and  (ii)  increase  £  by  1  from  £c(kl,l)  to  £n. 

For  each  value  of  the  pair  (£n,i)  calculate 
b(kl,£n,i)  using  equations  (2.92)  through  (2.97). 

For  each  value  of  (£,£n,i)  calculate  d(kl,£,£n,i) 
using  (2.91). 

Store  the  values  of 

(i)  d(kl,£,£n,i)  for  £n  <_  i  £  £c(kl,F), 

£c(kl,l)  £  n  _<  £c(kl,F)  and  £c(kl,l)  <  i  <  £n, 
and 

(ii)  dl(kl,£,i)  which  are  d(kl,£,£F,i) ,  (£F*»£c(kl,F)) , 
for  £c(kl,l)  ±  i  <_  m,  and  £c(kl,l)  <^  £  <_  £c(kl,F). 

If  kl  =  2,  go  to  Step  5.  Otherwise  go  to  Step  4. 
Increase  k'  in  increments  of  1  from  2  to  kl-1.  For 
each  value  of  k',  increase  £  by  1  from  £c(kl,l)  to 
£c(kl,F)  and  £'  by  J  from  £c(k',l)  to  £c(k',F).  For 
each  value  of  (k',£',£)  calculate  T(k' ,kl, £' ,£)  using 
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(2.90).  Store  all  values  of  T(k' ,kl,£* ,  Jl) . 

Step  5:  Increase  Jin  by  1  from  Jlc(kl,l)  to  Jlc(kl,F).  For 

each  value  of  Hn,  calculate  YCkl.Jln.Jln)  using  (2.138). 
For  each  value  of  Jin  >  Jlc(kl,l),  decrease  Jl  by  1  from 
Jln-1  to  Jlc(kl,l),  and  for  each  value  of  (H,Jln),  cal¬ 
culate  Y(kl,il,Jln)  using  (2.128). 

Store  the  values  of  Y(kl,Jl,JlF)  for  Jlc(kl,l)  <_ 

£  <  Jlc(kl,F)  where  JIF  -  Jlc(kl,F). 

Step  6:  Set  kl  “  kl+1.  If  kl  £  K,  go  to  Step  3.  Otherwise 
go  to  Step  7. 

Step  7:  Increase  k2  by  1  from  2  to  K.  For  each  value  of  k2, 
increase  Jl  by  1  from  Jlc(k2,l)  to  Jlc(k2,F).  For  each 
value  of  (k2,Jl),  calculate  H(k2,Jl)  using  (2.127). 

Store  all  values  of  H(k2,Jl). 

Step  8:  Increase  j  by  1  from  1  to  m  and  for  each  value  of  j 
calculate  Wj  using  (2.141).  Store  all  values  of  W^. 

Step  9:  Calculate  A  t ^  using  (2.140). 

Step  10:  Calculate  the  values  of  ,  j  -  l,2,...,m  using  (2.142). 

2.6.3  Algorithm  for  Calculating  E(B^)'s  and  p^'s 

The  values  of  the  E(B^)'s,  i  =  1,2,...,K,  are  calculated  using 
(2.15).  After  obtaining  the  value  of  E(I)  from  (2.10),  p^s, 
i  =  1,2,...,K,  are  calculated  using  (2.9).  These  calculations  are 
summarized  in  this  section. 
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For  computational  convenience,  equation  (2.15)  can  be  re¬ 


written  as 


E(B  )  =  1/ y [PROD(i)  E(I)  +  A  , SUM(i) ]  (2.143) 

l  nrt-i 


where 


PROD(i)  =  NX 


(2.144) 


E(I)  =  [  E  PROD(i)]"1  (2.145) 

i=l 

as  per  (2.10)  and  (2.144),  and 

Ni  K 

SUM(i)  -  [  E  (  E  {  E  A(v) }) ]  (2.146) 

v  =1  £=i+l  v  *0 
i  £ 

such  that  Vj  **  0  for  all  j  <  i. 

As  noted  in  section  2. 4. 2. 3,  A^| ^  corresponds  to  the  state 
(0,0,..., 0).  Among  the  quantities  required  to  calculate  E(B^) 
in  (2.143),  SUM(i)  is  the  only  one  which  involves  the  steady  state 
departure  probabilities.  SUM(i)  can  be  defined  as  the  probability 
that,  under  steady  state,  the  service  time  is  that  of  a  class  i 
customer  immediately  after  a  departure.  Because  of  the  way  in 
which  the  row  states  are  arranged  in  (1-P)\  SUM(l)  is  the  sura  of 
the  steady  state  departure  probabilities  corresponding  to  the 


first  N  {  n  (N  +1) } -1  row  states  of  (I-P)i,  SUM(2)  is  the  sum  of 
1  1=2  1  K 

those  probabilities  corresponding  to  the  next  N?{  II  (N.+l) }  row 

1-3 

states  and  so  on  with  SUM(K)  being  the  sum  of  the  steady  state 

departure  probabilities  corresponding  to  the  last  N  row  states 

K 

of  (I-P)1. 

Now  the  algorithm  can  be  summarized  in  the  following  steps. 

Step  1:  Increase  i  by  1  from  1  to  K  and  calculate  PROD(i)’s 
using  (2.144)  for  each  value  of  i.  Calculate  E(I) 
using  (2.145). 

K 

Step  2:  Set  i  =  1,  j  =  1,  and  j,  =  N  {  II  (N  +1)}-1.  Go 

1  1  1  1=2 

to  Step  3. 

J2 

Step  3:  Calculate  SUM(i)  =  I  A..  Set  i  =  i+1  and  j  -  j_+l. 

j=jx  3 

If  i  £  K,  go  to  Step  4.  If  i  >  K,  go  to  Step  5. 

K 

Step  4:  Set  j  „  =  j9  +  N.{  R  (N  +1) ) .  Go  to  Step  3. 

l=i+l 

Step  5:  For  i  =  1,2,...,K,  obtain  the  values  of  E(B^)  using 
(2.143)  and  then  the  values  of  using  (2.9). 

Step  6:  Stop. 

Using  the  values  of  the  P^'s  obtained,  the  mean  values  of 


the  system  performance  measures,  W  ,  L.,  W.,  L  ,  and  p  can  be 

qi  1  1  qi 

calculated  using  equations  (2.3)  through  (2.6). 


2.7  VERIFICATION  AND  COMPUTATIONAL  ASPECTS  OF  THE  ALGORITHMS 


2.7.1  Verification 


The  algorithms  developed  in  this  research  were  verified  using 

simulation.  The  algorithms  were  coded  in  Fortran  language  and  run 

on  IBM  370/155  system  to  obtain  the  required  output  measures.  A 

set  of  test  cases  were  chosen  with  different  values  of  the  input 

parameters,  K.,  N^,  vk,  X^,  (i  =  1,2,...,K),  and  exponential  service 

times  for  the  purpose  of  verification.  The  simulation  was  coded  in 

the  GPSS/360  language  and  run  on  the  IBM  370/155  system.  Point 

estimates  and  their  95%  confidence  intervals  were  obtained  for 

W  using  the  method  of  batch  means  [SARG  79]  with  5  batches.  The 
^i 

results  of  both  the  algorithms  and  simulation  for  different  test 
cases  are  given  in  Table  2.1.  It  can  be  seen  that  the  results 
agree  extremely  well,  thereby  verifying  the  algorithms.  For  the 
purpose  of  gaining  insight  into  the  computational  aspects  of  the 
algorithms,  some  other  test  cases  were  also  run  and  the  results 
are  given  in  Table  2.2.  These  are  discussed  in  the  next  subsec¬ 
tion. 


2.7.2  Computational  Aspects 

On  the  whole,  the  algorithms  behave  reasonably  well  and  are 
efficient.  This  is  especially  true  with  respect  to  those  algo¬ 
rithms  which  invert  the  (I-P)1  matrix  and  calculate  the  steady 
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CLASS 

i 

Ni 

A. 

l 

.  -  - 

W 

qi 

SIMULATION 

ALGORITHMS 

MEAN 

95%  C. 

I. 

i 

3 

1.0 

2.0 

0.65AA 

0.6A19 

to 

0.6669 

0.6593 

2 

3 

1.5 

2.5 

3.1695 

3.09A9 

to 

3.2A41 

3.1588 

3 

2 

3.0 

4.0 

23.6628 

21.2053 

to 

26.1203 

23.7393 

1 

5 

0.5 

12.0 

0.07A2 

0.0708 

to 

0.0777 

0.075 

2 

5 

0.6 

10.0 

0.1167 

0.1098 

to 

0.1236 

0.1176 

3 

5 

0.7 

11.0 

0.2226 

0.2056 

to 

0.2396 

0.2308 

1 

7 

0.1 

5.0 

0.0882 

0.0813 

to 

0.0951 

0.0895 

2 

7 

0 .  A 

10.0 

0.1297 

0.123A 

to 

0.1360 

0.1317 

3 

7 

0.5 

10.0 

0.2A90 

0.23A5 

to 

0.2635 

0.2592 

1 

15 

0.1 

2.0 

1.5052 

1.A198 

to 

1.5906 

1.5091 

2 

15 

_ 

0.06 

_ 

1.0 

22.2951 

20.A733 

to 

24.1169 

22.4878 

1 

1  j  0.5 

5.0 

0.1871 

0.1707 

to 

0.2035 

0.1976 

2 

2 

0.6 

7.2 

0.2371 

0.2197 

to 

0. 2545 

0.2418 

3 

2 

0.  A 

2. A 

0.2737 

0.2515 

to 

0.2959 

0.2794 

4 

2 

0.A5 

3.6 

0.A693 

0.AA68 

to 

0.4918 

0.4848 

5 

2 

0.30 

6.0 

0.7818 

0.7A31 

to 

0.8205 

J 

0.8134 

RESULTS  OF  VERIFICATION 

TABLE  2.1 
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state  departure  probabilities.  This  is  because  these  algorithms 
were  developed  taking  into  consideration  the  special  structure 
of  the  matrix.  To  illustrate  this,  the  number  of  additions  and 
subtractions  and  the  number  of  multiplications  and  divisions 
required  by  the  algorithms,  denoted  as  A1  and  A2  respectively, 
and  by  the  Gaussian  elimination  method  [HADL  61],  denoted  as  G1 
and  G2  respectively,  to  calculate  the  stationary  departure  prob¬ 
abilities  are  compared  in  Table  2.3  for  certain  cases.  It  can 
be  seen  that,  in  general,  the  number  of  operations  required  by 
the  algorithms  is  less  than  that  required  by  the  Gaussian  elimi¬ 
nation  method. 

The  number  of  operations  required  by  the  algorithms  depends 
on  the  number  of  classes,  the  number  of  states  of  (I-P)1,  and 
the  number  of  unity  elements  above  the  main  diagonal  of  (I-P)\ 
The  number  of  operations  required  by  the  algorithms  increases 
as  the  number  of  unity  elements  increases.  Therefore,  the  re¬ 
duction  in  the  number  of  operations  required  by  the  algorithms 
in  comparison  with  the  Gaussian  method  is  more  in  the  cases  where 
the  ratio  of  U,  the  number  of  unity  elements  above  the  main  di¬ 
agonal,  to  m,  the  total  number  of  states  in  the  matrix  (I-P)\ 
is  low  and  less  in  the  cases  where  this  ratio  is  high.  This  is 
illustrated  in  Figures  2.6  and  2.7  where  the  ratios  of  the  num¬ 
ber  of  additions  and  subtractions  and  the  number  of  multiplica¬ 
tions  and  divisions,  respectively,  required  by  the  algorithms 
to  the  corresponding  numbers  required  by  the  Gaussian  elimina- 
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RATIO  A2/G2 


0.8 


RATIO  U/m 


COMPARISON  OF  ALGORITHM  WITH  GAUSSIAN -I 
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tion  method,  that  is,  Al/Gl  and  A2/G2,  are  plotted  against  the 
ratio  of  U,  the  number  of  unity  elements  above  the  main  diagonal 
to  m,  the  number  of  states  of  (1-P)^  matrix  for  some  different 
values  of  m.  The  maximum  possible  value  of  the  ratio  U/m  is 
0.5.  Even  at  this  ratio,  the  reduction  in  the  number  of  opera¬ 
tions  required  by  the  algorithms,  as  compared  to  the  Gaussian 
elimination  method,  is  significant.  It  is  also  observed  that 
for  the  same  value  of  K,  the  number  of  classes,  and  the  same 
value  of  the  ratio  U/m,  the  reduction  increases  as  m  increases. 
Also,  in  addition  to  this  reduction  in  the  number  of  operations, 
the  storage  requirement  of  the  algorithms  is  always  less  than 
that  required  by  the  Gaussian  elimination  method,  as  the  whole 
(I-P)^"  matrix  need  not  be  stored  for  the  algorithms. 

When  the  number  of  states  becomes  large,  i.e.,  greater  than 
or  equal  to  250,  some  computational  difficulties  were  experienced 
with  respect  to  the  generation  of  the  elements  of  the  transition 
probability  matrix  P.  The  main  problem  is  due  to  the  floating¬ 
point  arithmetic  errors  in  computing  the  elements  of  P  as  per 
the  combinatorial  equations  obtained  in  Appendix  A.  This  prob¬ 
lem  was  minimized  by  using  double  precision  arithmetic  and  log¬ 
arithmic  operations  in  the  cases  run.  In  the  case  of  exponential 
and  hyper-exponential  service  time  distributions,  this  problem 
can  be  further  reduced  by  replacing  the  summation  of  the  largest 
range  with  a  single  product  term.  It  is  based  on  the  following 


combinatorial  relation  [KNUT  68] 


M  (-Ph-  „ _ l _  . 

:Qihl(h+v)  v[u+v 


(2.147) 


Another  problem  is  the  amount  of  computational  time  required 
to  calculate  the  elements  of  the  transition  probability  matrix,  P. 
Using  a  numerical  method  to  evaluate  the  integrals  necessary  in 
the  calculation  of  these  elements  may  reduce  such  problems.  Fur¬ 
ther  work  on  this  problem  is  desirable. 

When  the  state  space  becomes  large  the  storage  requirement 
also  increases,  especially  for  the  algorithms  which  invert  the 
(I-P)^  matrix  and  calculate  the  steady  state  departure  probabili¬ 
ties.  This  can  be  handled  by  using  direct  access  auxiliary  stor¬ 
age  devices  such  as  discs  or  drums.  Care  should  be  taken,  how¬ 
ever,  to  modify  the  program  such  that  the  input/output  time  does 
not  increase  significantly. 
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CHAPTER  3 


EXTENSIONS  AND  CONCLUSIONS 


3.1  INTRODUCTION 

In  this  chapter  we  first  discuss  the  extension  of  the  al¬ 
gorithms  for  obtaining  the  marginal  and  joint  time  average  prob¬ 
abilities  of  finding  a  certain  number  of  customers  of  different 
classes  at  the  facility.  Then  the  possibility  of  using  these 
algorithms  to  cover  mixed  class  models  with  infinite  capacity 
source  for  some  classes  and  finite  capacity  sources  for  the  other 
classes  is  discussed.  Finally  conclusions  and  suggestions  for 
future  research  are  given. 

3.2  MARGINAL  AND  JOINT  TIME  AVERAGE  PROBABILITIES 

This  extension  is  based  on  the  application  of  Level  Crossing 
Analysis  [SHAN  80]  to  a  special  case  of  discrete  state  processes 
known  as  piecewise  Markov  processes  [KUCZ  73] .  The  required 
equations  are  derived  and  then  modified  to  suit  implementations 
as  algorithms  in  the  following  subsections. 

3.2.1  Marginal  Time  Average  Probabilities 

The  aim  here  is  to  obtain  the  values  of  q^(n^),  for  i  ■  1,2, 
...,K,  and  n^  =  0,1,2,...,N  ,  which  are  the  marginal  time  average 


probabilities  of  finding  n^  customers  of  class  i  at  the  facility. 
In  order  to  understand  the  basis  of  the  approach  used  to  find 
these,  it  is  necessary  to  consider  the  sample  path  of  class  i 
customers  which  gives  the  number  of  class  i  customers  at  the 
facility  at  any  time.  It  is  shown  in  Figure  3.1. 


NUMBER  OF 
CLASS  i 
CUSTOMERS 


SAMPLE  PATH 


FIGURE  3.1 


An  arrival  of  a  class  i  customer  at  the  service  facility  is 
represented  by  an  upward  jump  and  a  departure  of  a  class  i  cus¬ 
tomer  from  the  service  facility  after  completion  of  service  by 
a  downward  jump.  A  level  n^  n^  <  N^,  is  considered  in  the  sample 
path  and  the  state  space  is  divided  into  2  mutually  exclusive  sets 
Set  1  consists  of  those  states,  less  than  or  equal  to  n^  and 
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set  2  consists  of  those  states  greater  than  or  equal  to  (n^+1) . 

Then,  as  per  the  Level  Crossing  Analysis  [SHAN  80],  r,(n. ), 

a  1 

which  denotes  the  rate  of  downcrossings  from  set  2  into  set  1, 
is  equal  to  r^n^),  the  rate  of  upcrossings  from  set  1  into  set 
2.  That  is. 


rd(ni>  “  ru(ni} 


(3.1) 


As  the  arrival  rate  of  class  i  customers  when  there  are  n^ 
number  of  customers  of  class  i  at  the  facility  is  equal  to 
^Ni_ni^i’  and  because  there  are  only  single  arrivals,  ru(n^)  is 
given  by 


ru(ni>  "  qi(ni)(Ni-ni)Xi  , 

for  i  =  1,2,...,K  and  n^  m  0,1,2, 


(3.2) 


-1  . 


Let  ni(n^)  represent  the  steady  state  probability  that  a 
departing  i  class  customer  leaves  behind  n^  customers  of  class 
i  at  the  service  facility,  (i  =  1,2,...,K;  n^  “  0,1,2,. .. ,Nj-l) . 
Then,  as  there  are  only  single  departures,  r^n^)  is  given  by 


rd(ni)  * 


ni(ni)piWi, 


for  i  *  1,2, ... ,K  and 


(3.3) 


n^  ■  0 , 1 , 2 , . . . , 1  , 
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where  gives  the  output  rate  of  class  i  customers.  From  (3.1) 

to  (3.3),  the  marginal  time  average  probability  of  finding  n^ 
customers  of  class  i  at  the  facility  can  be  obtained  as 


i  n  <ni)piyi 

q  (ni)  =  (Ni-ni)Xi  *  for  1  =  1»2 . K  and 


n^  “  0,1,2,..., N^-l  . 


q  (N  )  can  be  calculated  from 


(3.4) 


Nrl 

q1(Ni)  -  1  -  Z  qi(n±)  . 

V° 


(3.5) 


In  equation  (3.4),  H  (n^)  can  be  evaluated  using  certain 
elements  of  the  transition  probability  matrix  P  and  the  elements 
of  the  steady  state  probability  vector  A.  As  per  the  notation 

St 

introduced  earlier,  the  conditional  probability  that  the  (n+1) 
departure  leaves  behind  u^,^, . . .  ,Ujj  number  of  the  corresponding 
class  customers  at  the  service  facility,  given  that  the  n^  de¬ 
parture  leaves  behind  *  *  * ’-^K  number  of  the  corresponding 

class  customers,  is  represented  as  P[X^  ^  -  (u^,  i^, .  •  • ,  v^)  |  Xjj  * 

(j^» j£» • • • » Jr) J *  the  row  vectors  n  and  i  represent  (u^.u^, . . . .u^) 
and  (j,,j„,...,j  ),  respectively,  then  this  conditional  probability 

can  be  written  as  P[X  , ,  =  u  I  X  -  j].  These  are  the  elements  of 

— n+J.  —  — n 

P.  In  a  similar  way  the  steady  state  probability  that  a  departure 
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leaves  behind  **  **^k  numbers  of  the  corresponding  classes 

can  be  written  as  A(jJ.  Then  ni(ni)  is  given  by,  for  i  -  2,3, 

• • • ,K, 


n1^) 


Ni  N2 

I  z 

u^=0 


Ni-1  Ni+1 


ui-i=0  ui+r° 


NK  Ni  Ni+1 


1  It  E  E  £  P[X  °ulx  -j]A(1)}  (3.6) 

V°  Ji-M1+1-o  JK-o  -”*1 


+  ^[XnH_i=u|Xn-0]A(0>}] , 


where  =  0  if  l  <  i  and  ut  =  n^  for  nJL  -  0,1,2 . N  -1  , 


and  for  i  ■  1, 


n1^) 


n2  n3 


U2=0  u3*=0 


nk  nx  n2 

£  [{  £  l 

Y°  V1  ^2*° 


+  {PtX^^^-OlACO)}]  ,  (3.7) 


where  ^  for  ^  «=  0,1,2, ...  ,(^-1) 


In  (3.6)  and  (3.7),  ()  represents  the  row  vector  containing  all 
zero  elements,  i.e.,  (0,0,..., 0). 
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3. 2. 1.1  Implementation  of  the  Algorithm 

The  calculation  of  Il^n^)  as  per  equations  (3.5)  and  (3.6) 
is  the  main  task  in  the  calculation  of  the  marginal  time  average 
probabilities,  q1(n^).  Certain  elements  of  the  transition  prob¬ 
ability  matrix  P  and  the  elements  of  the  steady  state  probability 
vector  are  necessary  for  calculating  Il^(n^).  It  may  be  recalled 
from  section  2.6  that  only  those  elements  of  P  necessary  for  the 
calculation  of  the  steady  state  probabilities  are  stored  and  that 
the  steady  state  probabilities  are  calculated  after  the  generation 
of  the  elements  of  P  and  inversion  of  (I-P)^  matrix,  as  per  the 
previous  set  of  algorithms  developed.  Hence,  in  order  to  imple¬ 
ment  the  algorithm  to  calculate  q*(n^)  ,  it  is  necessary  to  store 
those  elements  of  P  necessary  to  calculate  Il^n^). 

In  equations  (3.6)  and  (3.7)  ,  the  vectors  u  and  j_  refer  to 
the  corresponding  column  and  row  states,  respectively.  Upon 
scrutiny  of  equation  (3.6),  it  can  be  seen  that  the  required 

elements  of  P,  which  are  the  conditional  probabilities  P[X  ,=u  I 

— n*T'i 

X  ]  correspond  to  the  row  states  which  has  zero  number  of  class 

1  customers.  Because  of  the  way  in  which  the  row  states  are 

arranged  as  per  the  steps  described  in  Section  2.4.2,  these  row 

K 

states  correspond  to  the  last  n  (N  +1)  rows  of  the  transition 

1*2 

probability  matrix  P.  Therefore,  it  is  necessary  only  to  store 

K 

the  elements  of  P  corresponding  to  these  last  H  (N  +1)  rows 

1*2 

as  they  are  generated  in  the  earlier  set  of  algorithms.  Equa- 
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tion  (3.7)  requires  the  elements  of  P  which  correspond  to  the 
row  states  having  the  number  of  class  1  customers  greater  than 
zero.  As  the  storage  of  these  elements  requires  additional 


space,  equation  (3.7)  is  to  be  modified  and  reduced  in  terms  of 

K 

only  the  elements  of  the  last  n  (N.+l)  rows  of  P.  To  achieve 

1=2 

this,  first  the  equilibrium  equation  (2.17)  is  written  as 


A(u)  = 


2 

I 


jr°  V0 


V° 


(3.8) 


where  I  u.,  I  j.  <  IN.  and  0  <  u.  <  N,  for  i  =  1,2 . K. 

i-1  1  i-1  1  i-1  1  -  i  -  i 

Equation  (3.8)  can  be  modified  by  rearranging  the  terms  and 
written  as 


1 

{  I 


N,  N, 


jr1  j^=o 


N„ 


.^[x^i-Hivi'iAa')  + 


P(Xu|1°u  |Xt=03A(Q) }  =  (A(u) 


(3.9) 


n2  n3 

I  I 

j2=0  33*° 


jK=0 


1  pix^uixrjJAan  • 


K  K  K 

where  I  j'  <  I  N  in  vector  j_’  and  j..=0  and  I  j  .?0  in  vector 

i=l  1  i*l  1  i=l  1 

X-  Using  (3.9),  (3.7)  can  be  rewritten  as 


n1(n1) 


n2  n3 

Z  I 

u2=0  u3=0 


nk 

Z  [A(u) 

Y° 


N2  n3 
z  z 
32=0  j3“° 


K 

Z 

V° 


(3.10) 


{P[4+1=u|4=jjA(i)}], 


3  =  0,1,2, ... ,N^-1  and  =  0  and  Z  j.  ^  0 


where  u,  =  n,  for 
in  j. 


K 

Z 

i=l 


Now  n  (n^)  for  i  =  1,2,...,K  and  n^  =  0,1,2 . N^-l  can 

be  found  by  using  (3.10)  and  (3.6)  in  terms  of  the  elements  of 
K 

the  last  n  (N  +1)  rows  of  P.  For  the  sake  of  computational 
11=2 

convenience  ,  equation  (3.10)  can  be  rewritten  after  interchanging 
summation  signs  as 


N„  N 


n1(n1) 


Sl(n1)  -  Z 


j2=0  j3=0 


NTJ 


Z  A(j^)Qi(n  ;j_)  ,  (3.11) 


jK=° 


where 


SK^) 


K 

Z  [ 


N, 


Z  A(u) ] 


^2  u£=0 


(3.12) 


and 


Q  (npi)  -  II  E  PIX^-ulV-Ul 


1=2  u^=0 


(3.13) 
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In  (3.11)  to  (3.13),  u1  -  for  n^  ■  0,1,2, .... (N^-l) ,  and  in 

1C 

(3.11)  and  (3.13)  j  =  0  and  Z  j.  >  0. 

1  £-2  * 

In  a  similar  way  equation  (3.6)  can  be  rewritten  as,  for 
i  -  2,3 . K, 


JI1^) 


{  Z 

V1  J 


N 


i+1 

Z 


i+1 


NK 

Z  A(j_)Q*(n.  jj.) } 

JR"0 


+  (A(0)Q1(ni;£) } 


(3.14) 


where 


Q  (nt;J)  -  i  l  z  Pix^^l^-L]) 

£=1  u  =0 

*i 

and 


g1^  ;o)  -  il  £  PQ^-ulx,-"))  • 

£“1  u  =0 

*i 


(3.15) 


(3.16) 


where  0  is  the  row  vector  containing  all  zero  elements.  In  (3.14) 

to  (3.16),  u^  =  n^  =  0,1,2 . N^-l  and  in  (3.14)  and  (3.15), 

j  >  1  and  j  ■  0  for  all  i  <  i  in  j_.  The  advantage  in  using 
equations  (3.11)  to  (3.16)  for  calculation  of  Hi(ni),  i  »  1,2, ...,K, 
is  that  the  storage  required  is  less  as  the  required  elements  of 
P  are  not  stored  individually  but,  rather,  as  sums  in  Q^Cn^^). 
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Now  the  algorithm  for  calculating  the  values  of  q^(n^)  can 
be  described  in  the  following  steps: 

Step  1:  Modify  Step  1  of  Part  A  of  Section  2.6.2  as 
follows: 

Calculate  and  store  Q*(n^;J[)  for  1  »  1,2,...,K 
using  (3.13)  and  (3.15)  ,  and  Q*(n^;£)  for  i  =  2,3,. ..,K 
using  (3.16). 

Step  2:  After  calculation  of  the  elements  of  A,  calculate 
Sl(n^)  for  n^  =  0,1,2, ... ,N^-1  using  (3.12). 

Step  3:  Calculate  the  values  of  Il^n^)  for  i  *  1,2,...,K  and 
n^  =  0,1,2, ... .Nj-l  using  (3.11)  and  (3.14). 

Step  4:  Calculate  q*(n^)  for  i  »  1,2,...,K  and  n^  ■  0,1,2, 

. . . ,N^  using  (3.4)  and  (3.5). 


3.2.2  Joint  Time  Average  Probabilities 

If  n  refers  to  the  row  vector  containing  the  elements 

(n^,n 

ing  ni>n2,,,,»IY  number  of  customers  of  the  corresponding  classes 
at  the  service  facility  can  be  represented  byq(n).  In  this  sec¬ 
tion,  expressions  for  obtaining  the  values  of  q(n)  for  i  *  1,2, 
...,K  and  n^  *  0,1,2,...,N^  are  derived. 

In  this  case  the  level  corresponding  to  the  state  n  *  (n. ,n., 
K  K  1  1 

,  E  n.  t  l  Nv,  is  considered  and  the  state  space  is 
i-1  1  i=l 

divided  into  two  mutually  exclusive  sets.  Set  1  consists  of  all 


2 . n^)  then  the  joint  time  average  probability  of  find- 
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i; ,  -I*-** 


states  n'  =  (n^,  n^ . n^) ,  such  that  0  £  £  ni  for  i  =  1,2, 

. . . ,K,  and  Set  2  consists  of  all  other  states.  Then,  because 
there  are  only  single  arrivals,  the  rate  of  upcrossings,  r^Cn), 
from  set  1  into  set  2  is  given  by 


r 

u 


(n) 


K  nl 

n2 

ni-i 

ni+l 

Z  A  .  (n . )  {  Z 

Z  ... 

z 

z 

i=l  1  1  n^=0 

n^=0 

n!  =0 
l-l 

ni+l=l 

where  n! 

i 

=  n .  in 
i 

n'  > 

°K 

£  q(n')} 

nK=0 


(3.17) 


and  X^Cn^)  =  (N^-nJA.^  for  i  =  1,2,...,K.  Also,  because  there 
are  only  single  departures,  the  rate  of  downcrossings ,  r^Cn),  from 
set  2  into  set  1  is  given  by 

K  . 

r,  (n)  =  Z  It*(n  )p.p.  .  (3.18) 

d  -  i=1  liii 

In  (3,18),  IT^n^)  is  the  steady  state  probability  that  a  departing 
class  i  customer  leaves  behind  n.  customers  of  class  i  and  less  than 

l 

or  equal  to  customers  of  class  1,  l?i,  at  the  service  facility  and 
gives  the  output  rate  of  class  i  customers.  Then  using  Level 
Crossing  Analysis  [SHAN  80]  as  r^n)  =  ru(n)> 

nl  n2  ni-l  ni+l  *¥. 

Z  Z...  Z  Z  ...  Zq  (ii ' )  } 

ni“°  n2=0  ni-i=0  ni+r°  ^=0  “ 


=  Z  Il*(n  )p  p 
i-1  1  1 


i  ’ 


K 

Z  A  (n  ){ 
i=l 


(3.19) 


Rearranging  the  terms,  equation  (3.19)  can  be  written  as 

K  K  . 

q(n)  £  A.n.  =  £  n1(n.)y.p.  - 

i=l  11  i=l  1  1  1  1 


K  nl  n2 

£  A . (n. ) {  £  £ 

i=l  1  1  n|=0  n2=0 


ni-l 

ni+l 

£ 

£ 

ni-l=° 

ni+l= 

K 

K 

£  n! 
3=1  J 

*  £  : 

j=l 

£  q(n' ) ) . (3.20) 

n^=0 


K 

In  (3.20)  the  elements  of  the  vector  n  are  such  that  £  n.  4 
K  j=l  J 

£  N..  q (N  ,N„ , . . . ,NV)  can  be  obtained  from 
j=l  J  1  1  K 


q(N1,N2, 


,Nk)  =  1  -  {  £  £  ...  £  q(n1,n2, . . . ,nK) }  . 

nl=0  n2=0  nK=0 

(3.21) 

K  K 

£  n.  t  £  N. 
j*ll  5-1  J 


Equation  (3.20)  gives  a  recursive  relation  for  finding  the 

s  t 

values  of  q(n)  starting  with  n  =  (0,0,..., 0),  which  is  the  (nri-1) 

row  state  of  (I-P)  matrix,  then  considering  the  row  and  so  on, 
s  t 

until  finally  the  1  row  state  given  by  n  =  (N^,N2> • • • ,  N^^jN^-l) 
is  considered.  The  values  of  Jl^n^  in  (3.20)  can  be  obtained 
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in  a  similar  way  as  described  in  Section  3.1.1.  The  algorithm 
for  finding  the  values  of  q  (n)  can  be  summarized  in  the  following 
steps : 

Step  1:  After  finding  and  storing  the  values  of  Il^Cn^) 
for  i  =  1,2,. ..,K  and  =  0,1,2, ...  ,IL-1,  set 
j  =  m+1 . 

Step  2:  Set  n  equal  to  the  jth  row  state  of  (I-P) .  Cal¬ 
culate  q(n)  using  (3.20). 

Step  3:  Set  j  =  j-1.  If  j  >  1,  go  to  Step  2.  Otherwise 
go  to  Step  4. 

Step  4:  Stop. 

3.3  MIXED  CLASS  MODELS 

In  this  section  the  possibility  of  using  the  algorithms  for 
mixed  class  models  with  infinite  capacity  sources  for  some  classes 
and  finite  capacity  sources  for  the  other  classes  is  considered. 
The  number  of  the  states  of  the  transition  probability  matrix 
becomes  infinity  if  one  or  more  of  the  sources  are  of  infinite 
capacity.  This  problem  of  infinite  number  of  states  can  be  e- 
liminated  by  assuming  maximum  limits  for  the  total  number  of 
customers  of  infinite  source  classes  at  the  service  facility. 
However,  care  should  be  taken  in  choosing  these  maximum  limits 
depending  upon  the  model  parameters  so  as  to  be  realistic. 

For  the  purpose  of  illustration,  a  mixed  class  model  with 
infinite  capacity  source  for  the  highest  priority  class,  i.e., 


class  1,  and  finite  capacity  sources  for  the  other  classes  is 
considered.  The  number  of  classes  is  equal  to  K.  The  capacity 
of  the  source  of  class  i,  i  =  2,3,...,K  is  assumed  to  be  N^, 
which  is  finite.  The  mean  time  interval  spent  by  a  class  i  cus¬ 
tomer,  i  =  1,2,...,K,  at  source  i  is  exponentially  distributed 
with  mean  1/A^.  Because  of  infinite  source,  the  probability  of 
a  certain  number  of  arrivals  of  class  1  customers  at  the  service 
facility  within  a  given  time  period  is  independent  of  the  number 
of  class  1  customers  already  present  at  the  facility.  This  is 
not  true  with  other  classes  for  which  this  probability  is  de¬ 
pendent  on  the  number  of  the  customers  of  the  corresponding 
classes  already  present  at  the  facility.  The  service  time  of 
class  i  customer,  i  «  1,2,...,K,  has  an  arbitrary  density  func¬ 
tion  with  mean  1/p^.  Class  i  customers  are  given  preference  over 
class  j  customers  for  service  if  i  <  j .  To  restrict  the  number 
of  states  of  the  transition  probability  matrix  from  being  in¬ 
finity,  it  is  assumed  that  the  maximum  number  of  class  1  cus¬ 
tomers  present  at  the  service  facility  at  any  time  is  M^.  In 
other  words,  there  will  be  no  arrivals  of  class  1  customers  in¬ 
to  the  service  facility  when  there  are  number  of  class  1 
customers  already  at  the  facility. 

This  model  can  be  analyzed  using  an  imbedded  Markov  chain. 

The  total  number  of  the  states  of  the  transition  probability 

K 

matrix  is  given  by  (M  +1) {  E  (N  +1) }-l.  For  example,  if  there 

i-2 
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10,  the  total  num- 


are  3  classes  with  =  20,  Nj  =  10,  and  = 
ber  of  states  is  equal  to  2540.  The  algorithms  for  ordering  the 
row  and  column  states,  inverting  the  corresponding  (I-P)^  matrix 
and  calculating  the  steady  state  departure  probabilities  are 
the  same  with  replacing  N^.  There  are,  however,  some  changes 
in  the  other  algorithms.  These  are  discussed  in  the  following 
subsections. 


3.3.1  Calculation  of  the  Elements  of  P 

The  probability  that  the  number  of  arrivals  of  class  i  cus¬ 
tomers  ,  i  *  2,3,...,K,  at  the  facility  within  a  time  period  t  is 
ra^,  when  there  are  n^  customers  of  class  i  at  the  facility,  is 
given  by 


p(mi:ni,t)  = 


N.-n. 
i  i 


-A. t  m.  -A.t  N.-n. -m, 
i.i,  l.iii 
i  )  (e  ) 


where  m . ,  n,  =  0,1,2,...,N.  and  m.  +  n.  <  N. .  This  is  the  same 
i  i  x  i  i  —  i 

as  expression  (2.20).  But  for  class  1  customers. 


p(m1;n1,t)  =  pGn^t) 

because  this  probability  is  independent  of  the  number  of  customers 
of  class  1  at  the  facility.  This  probability  is  given  by 
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pCn^Jt) 


e 


(3.22) 


-X  t  m 

axt) 


i 


Let  the  vector  j_ which  gives  the  state  of  the  process  at 
th 

the  n  departure  epoch  be  such  that  it  has  at  least  one  non¬ 
zero  element,  j^,  (l  =  1,2,..., K),  with  *  0  for  all  i  <  l, 
if  1.  In  other  words,  at  the  n1"*1  departure  epoch,  there  is 
at  least  one  customer  of  class  J l  present  at  the  facility  without 

any  customers  of  higher  priority  classes  waiting  for  service, 
s  t 

Then  the  (n+1)  service  time  is  that  of  a  class  l  customer. 

(i)  If  i  =  1,  then  the  (n+l)St  service  time  is  that 

of  a  class  1  customer.  If  the  vector  u  gives  the 

s  t 

state  of  the  process  at  the  (n+1)  departure 
epoch,  then  as  per  the  discussion  in  section  2.3, 


P[Xun»ulXu=j  1  =  o/“p(u1-j1+l;s1){  n  p(u1-ji;j1,s1)}dF(s1)  , 


(3.23) 


where  F(s^)  is 
service  time. 

p(ui"^i; Ji,sl) 
as 


the  distribution  function  of  class  1  customer’s 
Substituting  the  values  of  p(uj-jj+l;Sj)  and 
and  simplifying,  equation  (3.23)  can  be  written 


Pl4+r 


|X -II 


6io 


/  0[dF(8l)  , 


(3.24) 


where 
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and 


s; 


.  -V»,.  ,urJi, 

fe  (AjSj)  '■ 


K 

n  (l-e 
1=2 
n 


-A. 


Vi'ji,  "i 

)  (e 


.  Ni~ul 

*)  1  n 


(3.29) 


M  ~xt\uih+1  -  WW1. 

(l-e  )  (e  )  ] 


If  the  nth  departure  leaves  behind  an  empty  facility,  i.e., 
if  in  vector  ji  =  0  for  i  =  1,2,...,K,  then  as  per  the  equa¬ 
tion  (2.26)  in  section  (2.3),  the  conditional  probability  is 
given  by 

K 

P[X  .»u|x  =0]  =  Z  {Probfidle  period  is  terminated  by  a  class  i 
— n+i  —  — n 

customerlProMX^^uIX^jJ )  , 


where  0  =  (0,0,..., 0)  and  in  j_ ,  =  1  if  l  =  i  and  =  0  if 

if  i.  In  this  model  the  probability  that  an  idle  period  is 
started  by  a  class  l  customer  is  given  by 

v,“i  *  J.W1  ”hen  l>-2 

and 


K 

A  [A.  +  Z  N, 1 A  ] 
A  x  il=2  iX 


-1 


when  l  =  1. 
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Therefore  the  conditional  probability  is  given  by 


n^+1=^0]  -  Q1  +  Q2  , 


(3.30) 


where 


,-l 


Q1  =  + 


(3.31) 


and 


where  =  1  and  =  0  if  1  >  1  in  ^ 


K  K 

Q2  -  ViIXl  +  • 


(3.32) 


where  j  ■  1  if  £  ■  i  and  =  0  if  £  j*  i  in  j_. 


3. 3.2  Calculation  of  E(B^)'s  and  p^'s 


The  expected  length  of  the  idle  period  in  this  model  is  given 


by 


K 


,-l 


E(l)  =  [X  +  I  NX] 
1  i=2  1  1 


(3.33) 


Therefore  the  values  of  E(Bj)'s  are  given  by,  if  j  *  1, 
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E^)  *  l/u1IX1(A1  +  E  N1Ai)'X  + 


N1  K  Nil 


A(0)  E  (  E  {  E  A(v)  }) ] 
v^=l  J^2  v^=0 


and,  if  j  =  2,3, .. . ,K, 


(3.34) 


E(B^)  =  l/iiJ[NjAj(A1  +  E  NiA1)“x  + 


Nj  K  N£ 


A(0)  E  (  E  {  E  A(v) )) ]  , 
Vj-1  *-J+l  v£=0 


(3.35) 


where  v^  »  0  for  all  £1  <  j  in  v. 

The  values  of  p^'s,  i  *  1,2,...,K,  can  be  obtained  from  (2.9) 
with  E(I)  given  by  (3.33). 


3.3.3  Calculation  of  Mean  System  Performance  Measures 


The  value  of  p  can  be  obtained  from  (2.7).  The  equations 

(2.3)  and  (2.4)  giving  the  values  of  and  are  valid  only 

for  i  =  2,3,...,K.  For  i  =  2,3,...,K,  and  are  calculated 

first  using  (2.3)  and  (2.4)  and  then  W.  and  L  are  calculated 

1  qi 

using  the  expressions  (2.5)  and  (2.6). 

Class  1  behaves  like  a  submodel  with  limited  waiting  space 
of  capacity  M^.  The  steady  state  probability  that  a  departing 
class  1  customer  leaves  behind  customers  of  class  1  at  the 
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facility,  Ii(n^)  can  be  calculated  for  *  0,1,2 ,... ,M^-1  using 
the  equations  and  the  algorithm  developed  in  section  3.1.1. 

Applying  Level  Crossing  Analysis  [SHAN  80],  the  time  average 
marginal  probability  that  there  are  n^  customers  of  class  1  at 
the  facility,  q^(n^),  can  be  obtained  from  the  following  rela¬ 
tions  for  n^  =  0,1,2, ... ,M^: 

1  (n^) 

qx(ni)  =  - —  ,  for  n1  ‘  0,1,2 . M^l  and  (3.36) 


Mrl 

q1(M1)  =  1  -  Z  q1(n1>  . 
n^=0 


(3.37) 


Then  the  mean  number  of  customers  of  class  1  waiting  at  the 
facility,  L^,  is  given  by 

“i  i 

L  -  £  V  (»!>  • 

n^=0 


(3.38) 


The  mean  waiting  time  of  a  class  1  customer  at  the  facility,  W^,  is 
related  to  by,  [GROS  74], 


W1  “  VX1[1  "  q  (Mj.)! 


(3.39) 


Once  L,  and  W,  are  calculated,  L  and  W  can  be  obtained  from 
11  ql  ql 


L  =  L-  -  p . 

q.  1  1 


(3.40) 
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t 


and 

W  -  W.  -  1/u.  ,  (3.41) 

respectively. 

Thus  it  is  possible  to  use  the  algorithms  developed  in 
Chapter  2,  with  minor  modifications,  for  a  mixed  class  model 
with  infinite  source  capacity  for  class  1.  One  main  difficulty, 
however,  is  the  large  number  of  states.  This  increases  the 
storage  requirement  and  the  computational  time. 

3.4  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

The  main  contribution  of  this  research  iB  the  development  of 
a  solution  procedure  which  contains  a  set  of  numerical  algorithms 
for  finding  the  mean  values  of  system  performance  measures  for  a 
multiple  finite  source  queueing  model  with  fixed  priority  service 
discipline.  The  system  performance  measures  considered  are  the 
utilization  of  the  server,  the  mean  waiting  time  of  the  customers, 
the  mean  number  of  customers  waiting,  and  the  mean  number  of  cus¬ 
tomers  at  the  source  of  different  classes.  The  set  of  algorithms 
was  also  extended  to  obtain  marginal  and  joint  time  average  prob¬ 
abilities  of  finding  a  certain  number  of  customers  at  the  facility. 

The  possibility  of  using  the  algorithms  to  a  mixed  class  model 
was  also  discussed.  The  algorithms  were  verified  using  simulation. 
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Imbedded  Markov  chain  analysis  was  used  as  the  basis  of 
solution  development  of  the  algorithms.  The  system  was  imbedded 
at  departure  epochs  so  as  to  obtain  a  discrete  state  Markov 
chain.  The  structure  of  the  transition  probability  matrix  of 
this  Markov  chain  was  then  modified  by  suitable  arrangement  of 
its  row  and  column  states  so  that  numerical  techniques  could  be 
used  to  invert  the  final  matrix  and  to  obtain  recursive  solu¬ 
tions  for  the  steady  state  departure  probabilities.  The  mean 
values  of  the  system  performance  measures  were  related  to  these 
steady  state  departure  probabilities  using  the  regenerative 
method  and  Little's  formula.  The  time  average  marginal  and 
joint  probabilities  were  related  to  the  steady  state  departure 
probabilities  using  Level  Crossing  Analysis. 

The  set  of  algorithms  developed  in  this  research  provides 
an  alternative  to  simulation  for  studying  and  analyzing  this 
model.  Using  these  algorithms,  exact  values  of  mean  system  per¬ 
formance  measures  and  time  average  probabilities  can  be  obtained 
which  give  a  clear  insight  into  the  behavior  of  the  model.  How¬ 
ever,  there  are  a  number  of  areas  related  to  this  research  and 
the  model  that  need  to  be  investigated.  These  are  discussed  in 
the  remaining  part  of  this  section. 

As  mentioned  in  section  2.7.2,  some  computational  difficul¬ 
ties  were  experienced  while  generating  the  elements  of  the 
transition  probability  matrix  using  the  algorithms  when  the 
state  space  becomes  large.  Numerical  methods  can  be  used  to 


evaluate  the  Integrals  used  In  the  calculation  of  the  elements 


of  the  matrix.  Work  In  this  direction  Is  necessary  to  study 
the  effect  of  numerical  methods  on  the  accuracy  of  the  results 
and  the  computation  time  required  as  compared  with  the  combina¬ 
torial  methods  used  in  this  research. 

As  pointed  out  in  Chapter  2,  it  was  decided  not  to  use 
global  balance  equations  of  the  system  to  find  the  mean  values 
because  of  the  very  large  state  space  required  as  compared  with 
the  state  space  related  to  imbedded  Markov  chain  analysis  and 
its  limitation  to  service  time  distributions  with  rational 
Laplace  Transforms.  Exploration  of  global  balance  equations 
is  necessary  to  find  out  whether  any  special  structure  of  the 
matrix  related  to  the  equations  exists  which  may  lead  to  an 
efficient  solution  method  for  these  types  of  service  time  dis¬ 
tributions. 

The  use  of  the  solution  procedure  developed  in  Chapter  2, 
for  mixed  class  models  was  discussed  in  section  3.3.  Further 
investigation  is  recommended  to  determine  if  new  algorithms 
can  be  developed  for  mixed  class  models  which  are  more  efficient 
than  those  contained  in  the  solution  procedure  developed  in  th:?s 
research. 

The  behavior  of  the  model  investigated  in  this  research 
should  be  studied  to  find  the  effect  of  different  parameters 
such  as  the  number  of  classes,  number  of  customers  of  each  class. 
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service  and  arrival  rates,  and  service  time  distributions  on 
the  system  performance  measures.  Investigation  in  this  direc¬ 
tion  may  help  in  the  design  of  systems  which  can  be  described 
by  this  model. 

Only  the  mean  values  of  the  waiting  times  of  the  customers 
of  different  classes  are  obtained  in  this  research.  In  some 
cases,  the  mean  values  alone  are  not  sufficient  to  give  a  true 
picture  of  the  model  behavior.  An  important  measure  is  the 
variance.  Further  work  is  recommended  to  obtain  the  variances 
of  waiting  times  of  different  classes  of  customers. 

There  are  many  useful  results  available  for  multiple  class 
infinite  population  models,  but  not  applicable  when  the  popula¬ 
tions  are  finite.  One  such  case  is  the  rule  which  assigns 
priorities  to  the  different  classes  so  as  to  minimize  the  total 
delay  cost  based  on  mean  values  [CONW  67,  JAIS  68,  and  KLEI  76]. 
But  in  the  case  of  finite  population  models,  no  such  rule  is 
available  because  of  the  nonexistence  of  closed  form  solutions 
for  the  mean  values  of  different  performance  measures.  Future 
work  in  this  area  is  necessary.  Dynamic  priority  service  rules 
have  been  applied  to  multiple  infinite  source  models  and  useful 
results  have  been  obtained  [KLEI  76  and  WOOD  78],  In  the  case 
of  multiple  finite  source  models,  analytical  or  numerical  methods 
should  be  developed  for  finding  the  values  of  system  performance 
measures  with  dynamic  priority  service  rules. 
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Approximate  methods  for  solving  complicated  queueing  models 
are  gaining  increased  recognition  when  such  methods  can  lead  to 
operational  solutions  which  result  in  valid  decisions.  The 
main  advantage  in  using  approximate  methods  is  that  less  compu¬ 
tational  time  is  required  compared  to  exact  analytical  and  num¬ 
erical  methods  or  simulation.  There  are  many  ways  in  which  the 
performance  measures  in  this  model  can  be  obtained  approximately. 
Investigation  in  this  area  is  likely  to  yield  useful  results. 
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APPENDIX  A 


General : 
(i) 


(ii) 

(ill) 

(iv) 


(v) 


(vi) 


(vii) 


LIST  OF  SYMBOLS 


The  value  of  a  summation  is  considered  zero  if  the 

upper  limit  is  less  than  the  lower  limit.  For  example, 

12 

E  (•)  =  0  if  12  <  il. 
i=il 

The  value  of  a  product  is  considered  one  if  the  upper 

limit  is  less  than  the  lower  limit.  For  example, 

12 

n  (•)  -  1  if  12  <  u. 

£=«.l 

Matrices  are  denoted  by  boldface  letters  (B,  C,  P,  etc.). 
Column  or  row  vectors  are  represented  by  boldface  or 
lower  case  letters,  with  a  line  underneath  (A,  £, 
etc. ) . 

The  elements  of  a  matrix  are  represented  by  lower  case 
letters  with  the  row  and  column  numbers  as  subscripts. 

For  example,  the  element  of  matrix  C,  in  row  i  and 


The  i*"*1  element  of  a  column  or  row  vector  is  represented 
by  adding  i  as  the  subscript  to  the  letter  which  represents 
the  vector  without  the  line  underneath.  For  example, 
is  the  iC^*  element  of  the  vector  Z. 

The  inverse  of  a  matrix  B  is  written  as  B  \ 
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Specific: 

1A±: 


V 


n^n.): 


p 


i' 


Mean  inter-arrival  time  of  a  class  i  customer. 

Mean  service  rate  of  a  class  i  customer.  (■ 
reciprocal  of  the  mean  service  time  of  a  class  i 
customer. ) 

Steady  state  probability  that  a  departing  customer 
of  class  i  leaves  behind  n.  customers  of  class  i 
at  the  facility. 

Proportion  of  time  the  server  is  busy  with  class 
i  customers. 


p: 


Utilization  of  the  server. 


A: 


C: 


C 


p:i: 


E(B): 


Steady  state  departure  probability  vector  of  the 
imbedded  Markov  chain. 

Row  vector  containing  the  first  m  elements  of  A. 
element  of  A. 

Steady  state  probability  that  just  after  a  depar¬ 
ture  the  state  of  the  process  is  v  *  (v^.Vj* . . . ,vK> . 

Left  triangular  matrix  of  size  mxm  containing  the 
elements  located  along  and  below  the  main  diagonal 

of  (I-P)^  matrix. 

Matrix  which  contains  all  elements  of  C  and  all  the 

unity  elements  above  the  main  diagonal  of  (I-P)^ 
from  left  up  to  and  including  the  one  being  considered 

in  the  i^  iteration  of  phase  p. 

Expected  length  of  the  busy  period. 


E(B.): 


Expected  length  of  the  busy  period  during  which 
the  server  is  busy  with  class  i  customers. 


E(I):  Expected  length  of  the  idle  period. 

f(s.):  Probability  density  function  of  the  service  time  of 

a  class  i  customer. 


F(st): 


Distribution  function  of  the  service  time  of  a  class 
i  customer. 


G: 


Inverse  of  C. 
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;P:i. 


I: 

(I-P)1: 


K: 

Ac(k,n) : 
Ac(k,F): 


An: 

IF: 

Ar(i,j): 


m+1: 


p  •  i 

Inverse  of  C^‘  . 

Identity  matrix. 

Matrix  of  size  mxm  containing  all  elements  of  (I-P) 
except  the  last  row  and  the  first  column. 

Number  of  classes  (sources) . 

1  th 

Column  number  in  (I-P)  of  the  nu  column  of  the 
ordered  column  set  k. 

Column  number  in  (I-P)^  of  the  last  column  of  the 
ordered  column  set  k. 

£c(k,n) 

Ac(k,F) 

Row  number  in  (I-P)^  of  the  jth  row  of  the  ordered 
row  set  i. 

Mean  number  of  waiting  customers  of  class  i  at  the 
facility. 

Mean  number  of  waiting  customers  of  class  i  in  the 
queue . 

Total  number  of  states  of  P. 


N^:  Total  number  of  customers  of  class  i. 

p(m^;n  , t):  Probability  that  the  number  of  arrivals  of  class  i 

customers  at  the  facility  within  a  time  interval  t 
is  m^  when  there  are  n^  customers  of  class  i  at  the 

facility  at  the  beginning  of  the  time  interval. 


P: 

q1(ni) : 


q(n) : 


r(kl): 


Transition  probability  matrix  of  the  imbedded  Markov 
chain. 

Time  average  marginal  probability  that  there  are  n^ 
customers  of  class  i  at  the  facility. 

Time  average  joint  probability  that  there  are 
n  =  (n^,n 

corresponding  classes  at  the  facility. 

Total  number  of  row  states  in  the  ordered  row  set 
kl  of  (I-P). 


2,...,n^)  number  of  customers  of  the 
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r' (kl) : 


R(ic(ksn)) : 


R(Jln) : 


U: 


X  ,, : 
~n+l 


X  • 
^n+1" 


Z: 


Total  number  of  row  states  in  the  ordered  row  set 
kl  of  (I-P)1. 

Location  of  the  unity  element  above  the  main  diagonal 
of  (I-P)1  in  column  £c(k,n). 

R(Jtc(k,n) ) . 

Service  time  of  a  class  i  customer. 

Total  number  of  unity  elements  above  the  main  diagonal 
of  (I-P)1. 

Mean  waiting  time  of  a  class  i  customer  at  the  facility. 
Mean  waiting  time  of  a  class  i  customer  in  the  queue. 

State  of  the  process  at  the  n^  departure  epoch. 
ith  row  state  from  top. 

8 1 

State  of  the  process  at  the  (n+1)  departure  epoch, 
til 

vc  column  state  from  left. 

Row  vector  containing  the  negative  of  the  last  m  ele¬ 
ments  of  the  last  row  of  the  matrix  (I-P) • 
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APPENDIX  B 


EVALUATION  OF  THE  INTEGRAL  IN  THE  ELEMENTS  OF  P 


The  integral  Is  given  by 


o 


'VF<\> 


</V(st)dv 


where 


,  Kn  n  4isLVuiul  -A*vurU+1, 

{  n  (l-e  )  (e  )  ) (1-e  )  (e  ) 

i=l 


B.l  EXPONENTIAL  SERVICE  TIME  DISTRIBUTION 


f(V  “  V  *  *£  >  ° 


0  ,  otherwise 


K  -A.s.  u.-j.  -X.s.  N.-u.  -A..S,  u.-j 

c,  ,  .  c®,  n  ,,  i  L  i  i,  i  £.s  i  iWl  ^  !!.  i.  i  J 

/  0  f (s  )ds  =  f  {  H  (l-e  )  (e  )  }(l-e  ) 

o  5,  £  l  o  ,  , 

1=1 


(e  )  Vte  ds( 


"Vi 


Using  the  substitutions  y  =  e  and  a.  =  —  ,1-1,2, 

Vl 

,K,  the  final  form  of  the  integral  is  given  by 
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°  v<vdst  ■ {  *  I  ,  «-»  E0  u,  )  <-» 

&  X 


(B.X) 


JK  u  -j  ’ 

Z  l  K  (-1)  y  b}} 

hK=0  K 


for  0  <  u.  <_  N  ;  u.  -  j.  ^.0;  if  i  *  1,2, . . .  ,K,  and  0  <_  u.  <_  N.-l; 
i  i  X  1  /  ^  JtJC 

Ua  -  +  1  >  0. 

In  (B.l),  b}  is  given  by 


bJ  ■  i^<Wv1+V  +  1flAiwruirti»  + 11_1- 


(B.2) 


B. 2  HYPO- EXPONENTIAL  SERVICE  TIME  DISTRIBUTION 


(knP„)  s  « 

f(V  -  -—-(Trin - •  s*  >  0 


=  0  ,  otherwise. 


-k.ii„s„  A.. 

z  z  z  i 

Using  the  substitutions  y  =  e  and  a.  =  —  — —  ,  i  =  1,2 . K, 

JTH 


and  the  integral  relation 


/'*'xm[log(l/x)]ndx  =  - -  ■■■  Xi~  if  (m+1)  and  n  >  0, 

(nH-l)n  1 
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the  final  form  of  the  integral  is  given  by 


r°°n  ^  I*  ^'J1  /ul"Jl 

f  0of(so)ds  “  {  I  (-1)  Z  . 

°  1  1  1  h  =0  1  \  i  h  -0  1  hl  i 

£  1 


VjR  lu vrh)  \  2 

r  (  K  Kj  (-1)  *  bp 


(-1) 


(B.3) 


v° 


for  0<^u,  £.N;  u  -  j.  ^  0  if  i  =  1,2,. 
1111  ft  £ 

In  (B.3),  is  given  by 


,K,  and  0  £  u^  <_  N^-l; 


K  -k 

^’VW^V +  ^WW’ -  »  *• 


(B.4) 


If  is  taken  as  1,  which  makes  f (s  an  exponential  density 

2  1 
function,  b^  becomes  equal  to  b^. 


B.3  HYPER-EXPONENTIAL  SERVICE  TIME  DISTRIBUTION 


_UUSil  _V1£28£ 

f(V  -  PU  ^le  +  P*2  v£2e  •  S£  >  0 


0  ,  otherwise  , 


where  P£1  +  P£2  -  1. 
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This  is  a  two-stage  hyper-exponential  density  function.  It 
is  possible  to  have  more  than  two  stages  which  is  a  simple  exten¬ 


sion  of  the  two  stage  function. 

Writing  the  integral  as  a  sum  of  two  integrals  and  using  the 

~ullsl  Xi  ~v!2sl  Ai 

substitutions,  y1  =  e  ;  a...  *  -  ;  y5  *  e  and  a.»  •  -  , 

W£1  z  z  USL2 

the  integral  becomes 


W1 

/°0  f(s  )ds  =  {  E 
°  *  *  *  h  -0 


ln0-j0+l\  h'1'h 


SL 

h„ 


ve3k(vJk 

tyO  1 


(-1)  £ 
tyO 


t-D^  b3> 


urji'i  hi 

\ )  w> 


(B.5) 


for  0  £  u  £  N  ;  u.  -  j.  ^  0  if  i  =  1,2,...  ,K,  and  0  <  \i.  <  N  -1; 
1111  ^  £  x.  x, 

un  -  V1  -  °- 

3 

In  (B.5),  is  given  by 


b 


3 

l 


■  mt;rrfVW1+V  + 


E  X  (N  -u.+b.)>  +  1] 
i-1 

n 


(B.6) 


+  + 


l  X.(N.-u.+h  )}  +  l]"1  . 
i-1 
H 


It  P£l  is  taken  as  1  and  -  p^,  then  ^(s^)  becomes  an  exponen- 
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tial  density  function  with  mean  l/y^  as  in  B.l.  In  that  case, 

3  1 

becomes  equal  to  b^. 


i 
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APPENDIX  C 


EXPRESSIONS  OF  THE  ELEMENTS  OF  THE  INVERSE  G 


k+l:F 


Expressions  of  g  '  from  (2.87)  are  substituted  into  (2.98). 
*■»  J 


(i)  For  1  jS.  i  <  £c(2.1): 


I*or  J  = 


.....  k+l  £c(k+3-kl,F) 

g?  r  =  g,  ,  -  £  {  £  g.  j h (k+l , k+3-kl ,  £ ,  i)  } 

1,3  1,3  kl=2  £*£<:• (k+3-kl , 1)  *3 


£c(k+2,F) 

I  dl(k+2,£\i)[gs, 

£’*£c(k+2,l)  1 


k+l  £c(k+3-kl,F) 

I  {  Eg,  h(k+l,k+3-kl, £,£’)>] 

kl=>2  £*£c(k+3-kl,l)  ,3 


£c(k+2,F) 

-  E  g,»  ,dl(k+2,£  ,i) 

1,3  £’-£c(k+2,l)  ,3 


k+l  £c(k+3-kl,F) 

l  {  E  g  [h(k+l,k+3-kl,£  ,i) 

kl=2  £«£c(k+3-kl,l)  4,3 


£c(k+2,F) 

E  dl(k+2,£,,i)h(k+l,k+3-kl,£,£,n  (C.l) 

£’*£c(k+2,l) 


Consider 


£c(k+2,F) 

A1  =  h(k+l,k+3,£,i)  -  E  dl(k+2,£’ ,i)h(k+l,k+3-kl,£,£') 

£'«£c(k+2,l) 


Substituting  relations  (2.99) 


k+1  £c( k'.F) 

A1  -  D(k+3-kl,£,i)  -  E  {  E  D(k’ , ££,i)T(k+3-kl, 

k’=k+3-kl+l  U«lc(k’,l) 


£c(k+2,F) 

E  dl(k+2,£' ,i) [D(k+3-kl,£,£') 
£'-£c( k+2,1) 


k+1  £c(k',F) 

E  {  E  D(k\££,£')T(k+3-kl,k', £,££)}]  . 

k'=k+3-kl+l  ££*£ c (k ’ , 1) 


Consider 


k+1  £c(k’ ,F) 

A2  -  0(k+3-kl,£,£')  E  {  E  D(k', ££,£')  T(k+3-kl,k', 

k’*k+3-kl+l  IWc(kM) 


As  £'  >  £c(k+l,F),  because  of  (2.89), 


D(k+3-kl,  l,  £')  =  dl(k+3-kl,  l,  V) 

for  kl  ■  2,3, .. . ,k+l 


and 


' ,  £,  U)  } 


,££)}. 
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D(k',U,Sl')  =  dl(k'  ,Zl,l')  for  k'  -  k+3-kl+l , . . . , 


k+1 . 


Therefore, 


k+1  £c(k',F) 

A2  ■  dl(k+3-kl,  *.,&')  -  Z  {  2  dl(k'  ,U,SL')  T(k+3-kl,*,M)}  . 

k '  =k+3-kl+l  zz=lc(k'  ,1) 


It  can  be  seen  that  the  expression  for  A2  is  identical  to  the  ex¬ 
pression  for  T (k+3-kl , k5 , £ , & 1 )  given  by  (2.90)  for  k5  =  k+2.  Hence 

A2  *  T(k+3-kl,k+2,J^.,Jl,)  . 


Now  A1  becomes 


k+1  i,c(k',F) 

Al  -  D(k+3-kl,*,i)  I  {  Z  D(k’,M,i)T(k+3-kl,k',£,M)} 

k'-k+3-kl+l  Jt«t,c(k',l) 


£c(k+2,F) 

Z  dl (k+2 , l ?,i) T (k+3-kl ,k+2 ,£,&')  . 

Z’mZc(k+2,l) 


Now  i  <  Zc (k+2, F)  and  so 


dl(k+2,*',i)  -  D(k+2,£',i) 


because  of  (2.89). 


Therefore, 


k+2  Hc(k',F) 

A1  **  D(k+3-kl,H,i)  Z  {  Z  D(k'  ,M,i)T(k+3-kl,k'  ,H 

k*k+3-kl+l  U-£c(k',l) 

Al  is  identical  to  the  expression  for  h(k7,k+3-kl,H,i)  given  by 
(2.88)  if  k7  is  taken  as  (k+2). 

Therefore, 


Al  =  h(k+2,k+3-kl,H,i)  . 


Now  equation  (C.l)  can  be  rewritten  as 


k+l:F 

gi,j  =  gi,j 


He (k+2, F) 

Z  g  ,  dl(k+2,H',i) 

H'«Hc(k+2,l) 


k+1  Hc(k+3-kl,F) 

Z  {  Z  g  h(k+2,k+3-kl,H,±) } . 

kl»2  H=Hc(k+3-kl,l) 


Since  i  <  Hc(k+2,F), 


dl(k+2,H',i)  -  D(k+2,H\i) 


as  per  (2.89). 

D(k+2,H',i)  is  identical  to  the  relation  for  h(k7,k+2, £' ,i) 

k+1  •  F 

as  per  (2.88)  when  k7  ~  k+2.  Therefore,  g,  ’  can  be  written  as 


,U)} 


...  -  k+2  Hc(k+4-kl,F) 

g*  ,  ■  g,  ,  -  I  i  £  gp  ,h(k+2,k+4-kl,H,i)} . 

1,3  1,3  kl-2  H=Hc(k+4-kl,l)  *3 


(ii)  For  He (2,1)  £i  <  Hc(k+1,F): 


For  j  “  1,2, ... ,m, 


k+2  He (k+3-kl, F) 

fc+l.i  _  _  E  {  j  g  h(k+l,k+3-kl,H,i) } 

1,3  kl=2  H=Hc(k+3-kl,l) 


He (k+2, F) 

E  dl(k+2,H',i)[g  , 

H '=>Hc(k+2, 1)  *  ’ 3 


k+1  Hc(k+3-kl,F) 

-  I  {  I  g  h(k+l,k+3-kl,H,H>)}]  . 

kl=2  H=Hc(k+3-kl,l)  ,3 


This  is  the  same  expression  obtained  in  (i)  for  1  <_  i  <  He (2,1) 


without  the  term  g.  ..  Therefore,  following  the  simplifications 
**  J 

k+l:F 

in  (i),  the  final  expression  for  g.  .  can  be  written  as 

i » j 


k+2  Hc(k+4-kl,F) 

=  -  I  {  Z  g  .h (k+2 ,kl+4-kl , H , i) }  . 

kl=2  H=Hc(k+4-kl,l) 
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